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1. Introduction. If a fluid is streaming past a body of revolution with a velocity at 
infinity parallel to the generators of the body, then it is well-known [2] that the boundary 
layer induced on it is equivalent to the boundary layer on a flat plate provided only 
that the thickness of the layer is small compared with the radius of the body. Recently 
an increasing amount of attention has been paid to the flow of compressible fluids over 
slender bodies of revolution at high speeds. Since the displacement thickness 6, of the 
boundary layer on a flat plate is of the order M’xR~’”’ where z is the distance from the 
leading edge, VM is the Mach number of the flow in the main stream and R is the Reyn- 
olds number based on «x, there is a real possibility that 6, may be at least of the same 
order as the radius of the body. It is of interest to examine the modification to the 
boundary layer which is then necessary. We restrict attention to incompressible flow, 
but it is hoped that the methods developed here may be extended to compressible 
boundary layers. 

The boundary layer on the outside of a circular cylinder near to the leading edge 
has already been considered by R. A. Seban and R. Bond [3]. They assumed that the 
axis of cylindrical tube was parallel to the direction of flow and that the boundary 
layer was of zero thickness at the leading edge of the cylinder. They expanded the 
stream function in a power series in ascending powers of (x/a)'”” where z is the axial 
distance from the leading edge and a the radius of the cylinder. The leading term is 
the Blasius solution for a flat plate in a uniform stream and they computed numerically 
the two following terms. Among other things they found that the skin friction coefficient 
C, on the cylinder is 


edt.) : (18) . (100) | 
0.664(-*-) E + 0.53(773) — 0.071\773) + -+- |, 


where v is the coefficient of kinematic viscosity and U the (constant) velocity in the 
main stream. The initial effect of curvature is therefore to increase the skin friction. 
It is the aim of this paper to examine the boundary layer at large values of x. For 
this purpose we are not concerned with the precise form of the body near z = 0, although 
in fact we shall assume it to be a circular cylinder for all x > 0, and little modification 
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would be needed to account for the effect of a rounded nose. We show that if the velocity 
of the main stream is constant, the skin friction may be expanded as a double power 
series in descending powers of log ¢ and £ of which the leading term is proportional to 
[log ¢]”" where £ = 22zv/Ua’. It is thus of larger order than at corresponding points of a 
flat plate, so that the increase found by Seban and Bond [3] near the leading edge is 
maintained a long way downstream. The thickness of the boundary layer is, however, 
only slightly reduced in comparison with the flat plate and this implies that most of 
the changes in velocity occur relatively close to the cylinder. 

The asymptotic form of the boundary layer is also considered when the main stream 
velocity is proportional to x” where m < 1. If m > 1 the boundary layer thickness 
diminishes as x increases so that the effect of curvature is small when x/a is large, and 
if m = 1 the boundary layer is of constant thickness and there is a solution dependent 
on one parameter only [1]. We show that if —} < m < 1 the leading term in the asymp- 
totic expansion of the skin friction is also proportional to [log ¢]~’ and we infer that 
the effect of the curvature of the body is to delay separation. 

The contrast between the boundary layers on a flat plate and on a circular cylinder 
is not confined to the skin friction, for on the flat plate the boundary layer, while diffusing 
outwards, retains a similar form, whereas on the circular cylinder the form is always 
altering and ultimately tends to that given by Oseen’s approximation. The difference 
between plane and axially symmetric boundary layers may also be illustrated as follows. 
Suppose we have a viscous fluid occupying the region between two parallel planes a 
distance d apart, of which one is at rest and the other moves with velocity U parallel to 
itself. Then the velocity u of the fluid is given by 


u = U(1 — y/d), 


where y measures distance from the moving plane. It is seen that the form of the velocity 
profile is independent of d, being always a straight line. Now consider the corresponding 
problem in axial flow. We have a pipe bounded by two concentric cylinders of which the 
outer of radius b moves parallel to itself with uniform velocity U and the inner of radius 
a is at rest. Then the velocity u of the fluid in the pipe is 

, log (r/a) 


“= wee 
log (b/a) 


where r measures distance from the axis. It is seen that the radii of the cylinders are of 
great importance in determining the form of the velocity profile. For example if we fix 
b and r, then as a — 0, u — U so that for small a there is a boundary layer in the vicinity 
of the inner cylinder. If we interpret these two solutions as the asymptotic form of the 
boundary layers when the fixed surface is semi-infinite, it may be deduced that in the 
two-dimensional problem the boundary layer on the fixed surface grows until it has 
spread uniformly throughout the fluid, while in the axially symmetric case with a small, 
it is always confined to the neighbourhood of r = a. 

2. The equations of the boundary layer. We consider a circular cylinder of radius a 
whose axis occupies the positive half of the z-axis. Let r denote distance from this axis 
and u, v the components of the fluid velocity along the directions of x and r respectively. 
Then, if v is the coefficient of kinematic viscosity, p the pressure, and p the density, the 
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equations of motion in the incompressible boundary layer are 


+ 








fF) re) 
Zz (ru) + or (rv) = 0, 
Ou ou ldp , va ( ou) 
— poe a Ete (rE a | 
a a= ode ret al (2.1) 
and o(1) = 12. 
p or 


The boundary conditions are 
u=v=0 when r = a,x > 0, 
u— U(x) asr— @ for fixed x and as x 0 + for fixed r > a, 

where U(x) is the velocity of the fluid in the main stream. The third equation may be 
interpreted as implying that the pressure variation is small across the boundary layer. 
We shall assume this to begin with, and then, using the solution we obtain, it may be 
shown that the pressure variation across the boundary layer is O(pUv/x) and negligible 
when « is large. The values of U in which we shall be interested in this paper will be 


proportional to x” so that 


_lop _ mU" (2.2 


p Ox = 
and from the equation of continuity we can define a stream function 


viV(E, 7), 


where 
§ = 2rv/Ua’ and n =r U/2x». (2.3) 

Then 

a 
u=12orm -u= (2.4) 
r or On 
and 
1 a , v ov ov 

= ——< av) = —] ¥ - (1 - mn — mg = 2.5 
t a (va) : E (J m)n an + (1 m)é ae |. (2.5) 


The equation satisfied by ¥ is 


pith. wate | . (22)'] me a (2 ov _ oF a4 . 
2n an +(2+ VW) a + m| 1 an (1 m)é an an ak aE onl” (2.6) 





The boundary conditions are 


ay ov ad 
—=y + (1 — mt =0 at nt = 1, (2.7) 
On dé 

av/dn — 1 as "> @ for allé and as t— 0+, 


provided né > 1. 
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We note that if m = 1 then 7 is a function of r only and that a solution may be found 
with WV a function of 7 only. This special case has already been discussed by Cooper and 
Tulin [1]. If m > 1, then § +0 as xz -—> ~. The boundary layer decreases in thickness as 
x —o so that the effects of curvature are greatest at the leading edge and may be neg- 
lected in the asymptotic expansion, which may therefore be reduced to that for a flat 
plate with the same mainstream velocity. 

3. Uniform mainstream. Here m = 0 and to obtain the first and crucial term of the 
asymptotic series for Y we assume that as § > o~, tdW/dt — O for fixed n. The right 
hand side of (2.6) is then zero in the limit  - ~ and the boundary conditions reduce to 


ov aw 
v= hes 0 at 9 = 0 and —— | as n> ©. 
on On 


Equation 2.6 may be integrated once to give 


ay ; & an | 
1 oy = A exp | J, 2n | (3.1) 


where A is a constant. We deduce that either d°V/dn? ~ 1/n near » = 0 which implies 
that the boundary condition on dV/dn at 7 = 0 is violated or A = 0 and @°W/dn’ 
except possibly at » = 0. We accept the second possibility and obtain as our first approxi- 
mation 


Vv = 7». (3.2) 
This is simply the stream function of the undisturbed stream and satisfies (2.6). The 
boundary condition at » = £‘ is invalidated however. The approximation may be 


improved by substituting back into (3.1) and applying the boundary condition on 
OV/dn at » = &* instead of at 7» = 0. We find that 








av ae Ae? 
7 an’ d ’ 
whence 
~ = A ee 
Now when 7 is small 
et dz _ (<P (2) 
I a F —to 08 > = + : (n + 1)(n + 1)! ‘ aie 
where log C = 0.577 --- and is Euler’s constant. 
Hence, since dV/dn = 0 at né = 1, 
= [log (2/C)]"* + O[€ log (2&/C)}", (3.4) 
the second approximation to dV/dn is 
1 a ie , 
1 = peerre / =, (3.5) 


and we note that the correction term just found is small when ¢ is large except near nf = 
This suggests that we write 


‘s > Fn) _ (=) 
vant 2 Tog @/OF + Neioge — 








1955] BOUNDARY LAYER IN AXIAL INCOMPRESSIBLE FLOW 117 


and investigate the properties of F,(n). If our expansion is to be meaningful then it must 
be possible to make the error in W as small as we please by taking ¢ large enough and 
taking a sufficient number of terms. Certainly this is not true if we retain the first term 
only, but we shall show below that the successive terms are of uniformly decreasing order 
which suggests that our requirement is met. 

It is of interest to note that the first two terms of this expansion are the same as 
those obtained on Oseen’s approximation when z is large and it is likely that a similar 
result would be found whatever the cross section of the cylinder. In contradistinction 
to the boundary layer on the flat plate, it appears in fact that Oseen’s approximation 
becomes more accurate as x increases. 

Substituting in the differential equation and comparing coefficients of (log 2¢/C)~* 
we obtain 
2nFi"’ + (24+ FY = -—(s — IF, — F,-,FI’ 

ont (3.7) 
— 2D (PP + PP — Fie], 
t=1 


where the primes denote differentiation with respect to 7. In the determination of the 
F, we shall take as boundary conditions 
Y=0 at »=0, dv/dn=0 at » =€' and dV/dn-> 1 as nom. (3.8) 


The first of these conditions is inaccurate and so later on we must investigate how large 
is the error thereby incurred. We now prove that if 


Fi(n) = D, log nC + E, + A:n(log $nC)’ + Bin log 4nC + O(n) (3.9) 


near 7» = 0, for 1 < ¢ < s — 1, where the A’s, B’s, D’s and £’s are constants, then 
F‘(n) has a similar form near 7» = 0. It may then be inferred that (3.9) is the correct 
form for F/(m) for all s > 1 since it is the correct form when s = 1. From (3.9) it follows 
that near 7 = 0 


F, = Din log $nC + (E. — D.)n + 4Acn*(log nC)’ + 3(B. — A.)n* log $nC + O(n’), 
since F, = 0 at » = 0, and 
Fi’ = Dn + A,(log $nC)* + (2A, + B,) log $nC + O(1). 
Hence substituting into (3.7) we find that 
QnF i" + (2 + n)Fi’ = 2A,(log $C)? + 2(B, + 4A,) log nC + O(1), (3.10) 


where 
s-2 
2A, = — >, tD,D,-1-1 
t=1 
and 


e-2 
—2B, — 8A, =(s — 1)D,-, + D,D,-, + ps (tD,D,-.-1 baa D,D,-.) 


+ FWD Ey aar + Di-s-sE)- 
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Integrating (3.10) once, we obtain 
nF!’ = A,n(log 4nC)° + (B, + 2A,)n log $nC + D, + O(n), 
whence 
F’ = D, log 4nC + E, + A,n(log 4nC)” + Bn log 3nC + O(n), 
and, comparing with (3.9), we see that the desired result is proved. Thus, near 7 = 0 


ov — > D, log An + E 


—! + Olé "flog (2&/C)]~’). 


On = {log (2&/C)| 
In particular, when 7 = 1/, dV/dn = 0, so that 


- D, E, . 
b= 2.3 PUCTI/6)) Ga Zo3 ay (3.11) 


(Ie 
yr it = 
Zé | 1 [log (2 


and therefore 


D, = 1: D, = E,-., >t. (3.12) 
We may now express A, , B, in terms of D, obtaining 
2A, = — > tD.D,-1-1 (3.13) 
and 
] 


B,+4A, = —(s— 1)D,_, — = D {2(t — I)D.D,_, + tD,Dy-1-1}. (3.14) 


Mz 


Since EF, is determined from the condition that F/ — 0 as n > o, and A,, B,, D, are 

dependent only on the E’s we may determine as many of them as we please by successive 

substitution. The procedure is illustrated below by the determination of D, and D, . 
The equation of F, is 


QnF 1’ + (2 + FY’ = 0, 


with solution 


ree Dee, red, [f ore 
and | : (3.15) 

Pie Dig fo" ~apa - «”. 
From (3.3) and (3.13) it now follows that 

dD, = and E, = 0. (3.16) 
The equation for F, is 
2nF 3” + (2 + oF?’ = —Fi — FFY’, 
whence 
© (ane Fy) = —1 F, — Fre” 
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so that 


One Fi! = -[ Fie”? dn — [ F,& 
, . (3.17) 
n/ ~f oe 
= (2? +542) [ « : = + 4 log 1 - Qe" + C2, 
“9 
where C, is a constant. 
When 7 is small the right hand side is equal to 
—4 log $nC + 4 logn —2+C, = 2D,. (3.18) 
Hence from (3.12) and (3.16) 
C, = 2 + 4 log $C 
and, since F; — 0 as 7 > &, 
; si ee oe ee 
E, = -3 | é = \(2e ++2) | € — + 4 log » — 2 24 C, 
lies ” =“ ‘ (3.19) 
= —2 log2 — = = D,. 


Further terms of the series may be found in a similar way although the computations 
soon become very complicated and recourse to numerical methods is then necessary. 

In order to complete this discussion of the asymptotic boundary layer on the circular 
cylinder we must now investigate the leading error term in (3.6) to show that it is 
sufficiently small for large x. We write 


‘i — ' 
Y=1+ DL in@ort? (3.20) 


and suppose that the F, are all known. The differential equation (2.6) is satisfied if 
@ = 0 but not the boundary conditions and, neglecting terms involving 7° and powers 
of log 4Cn, it follows from (3.9) that near 7 = 0 
Fi = D, log nC + EB, + A,nllog 1nC)? + B.n log $nC + Con +::: 
where the A’s, B’s and C’s may be expressed in terms of the D’s and E’s. Hence when 
n = &' we have, neglecting ¢°’, 
A, . B, C. 


o®P 
dn ~ 2 iiog Gt anal 2. Flog 2 (2¢/C)]"* 2 iflog (2¢/C)]" 








and 

+t Lae awort Llegaor | O20 
Now from (3.13) and (3.14) and the solution for F, already worked out, 

A,=0, A,=0, A,=-3, A,=0 

B, = 0, B, = —1, B, = 0, 

C= 3 Gut 
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and hence the boundary conditions for ® are 


ae 1 5 


dn [log (2/C)] flog (2E/C) fe vee, 





o® Oo® 2 2 
= = _s __* = ‘ wie 22 
(e} ” on E[log (2&/C)] E[log (2¢/C)]}" + (3.22) 


+E 


when 7 = £' and £ is large, and 06/dn7 ~ 0 as n-— ~. 
If now we put 


G,(n) ( 1 ) 
6 = —— _4 + oi ——-—_ 
E[log (2¢/O)F * “\éflog (28/0) 57” 
then 27G{’’ + 2G{’ + 7Gi’ = 0 
with the solution 
v re 
Gi(n) = | "aia = (3.23) 


Jn 


Further terms of the series for may be obtained if desired but sufficient has been 
done to show the form which the asymptotic expansion of d¥/dn must take, viz., 


ov ~ P? ils) : 
Ka1+ 5 _, (3.24) 
an 2s [log (2E/C)] 

where the P’s are to be determined successively first for s = 1 and all ¢, then s = 2 and 


all ¢, etc. It is noted that there are no terms in the expansion of dV/dn with ¢ = 0 owing 
to the form of the expansion in (3.3). Further it is only necessary to obtain the P,,,(7) 
for 1 < ¢ < p in order to determine P,,,-. although such terms are of higher order 
than any of the P,,, . Finally it is pointed out that as soon as we begin to examine the 
P,., the effect of the shape of the body near x = 0 becomes of importance so that for 
example a finite shift of the leading cross section of the cylinder along the z axis will 
change them. 
The skin friction on the cylinder 


au) il (ov 
(24) . = pt x a). e - 
é 3.20) 
9 a [ F{’(n) F3'(n) ; , G1’(n) 
= pU* - (De /™ eT Ps a/v + °° |: 
x Llog (2&/C) ° [log (2€/C)] é[log (2&/C)] 


Now when 7 is small 


Fi'(n) = 1! — 4 + O(n); F2'(n) = —log nC/2 + § + O(n); 
Fi/(n) = —n ‘(2 log 2 + 4°) — 4(log 4nC)* — log $nC + O(1); 


Fi/(n) = 1 'Dz + O(log n), and Gi’(n) = n* + O(1). 


Hence after some reduction 


Ou 2uU | 2 log 2 + 3r° 
(. ou) a E Arv/Ca?U) (log (has CaO) P ieee 
= (3.26) 


a’ U Z 
+ Qxv ‘3 log (4av/Ca,U) me 7 + w. 
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For a flat plate at a distance x from the leading edge the skin friction is 
[u(du/dy)],-0 = 0.332nU(U/va)"”’, 


which is of smaller order than (3.26) when zv/a’U is large. Thus the trend noticed by 
Seban and Bond [3] near the leading edge of the cylinder is emphasized at large distances 
downstream and we may conclude that the effect of the curvature of the body is to 
increase the skin friction, especially at large distances downstream when the boundary 
layer is thicker than the radius of the cylinder. The skin friction is then almost constant, 
diminishing like (log x)~’. 

We may also calculate the velocity profile at large distances from the cylinder when 
x is large and we find that 


— 
u xi 


1 y 
U-~ ~~ 9 log 2&/C” 


The boundary layer thickness is ultimately of order 





1/2 
v 
“Gr log +. ans) : 


so that the effect of curvature on it is small thinning it out by a factor (log x)~'”’. 


4. Main stream velocity proportional to x”. There is a difficulty with the boundary 
layers caused by main stream velocities of this kind at the leading edge of the cylinder, 
especially if m < 0. We are however interested in the behaviour of the boundary layers 
at large distances downstream, so that we may avoid the difficulty by supposing that 
the main stream approaches this form at large distances downstream and our solution 
will give the asymptotic form for the ensuing boundary layer. 

As pointed out in Sec. 2 our method is only applicable if m < 1 and then we write 


= F,. 
Want one, (4.1) 
«t=1 (log &) € 
with the boundary conditions 


ov ov 
~~" as 7 ©, ag UT me 


aw 


ae = 0, 


when 7 = &"'. 
The values of P,,, may be obtained seriatim. 

Of particular interest is the effect of curvature on the tendency of the boundary 
layer to separate, and some information on this score may be obtained from the leading 
terms of (4.1). If we set P,.,(n) = H(n) then 


2nH’"” + (2 + »)H” — 2mH’ = 0, 
with boundary conditions 
H’ = —logét at n=" and H’—-0 as n> @, 


The appropriate solution is 


, * «alt= 1y" 
H ()= A | ‘ n(t=1) " (4.2) 
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if 2m > —1. Defined by this integral, H’(n) — 0 exponentially as 7» © and near 7 = 0 


ol 


H'(n) = —A log} nC + A | a (1 — o)'"—1)+- 


0 


Hence A = —1 and 


9TT 
(2) = ae + higher order terms. 
If 2m = —1 the integral is singular and in fact the method breaks down because H’ = 
Ae~”’ and does not behave logarithmically near 7 = 0. It is not clear at the moment, 
what is the correct procedure in this case. Provided —2m is not a positive integer, the 
method outlined in the previous section is formally adequate but in view of the break- 
down at m = —# we cannot assume at present that it is physically significant. 

Now if the main stream along a flat plate is proportional to x” we know that there 
are similar solutions of the boundary layer equations with positive skin friction if 
m > —0.0904 while from our work above we see that on a circular cylinder the skin 
friction is certainly ultimately positive if m > —0.5. We must bear in mind however 
that the conventional boundary layer comes to an end at separation, and that the 
same may well be true of a circular cylinder. There is a real danger that separation will 
occur at a finite value of x if —0.5 < m < —0.0904. Nevertheless we may conclude 
that the effect of the curvature of the cylinder, when the boundary layer has a thickness 
comparable with its radius of curvature, is to delay separation. 
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SOLUTION OF LINEAR EQUATIONS BY DIAGONALIZATION OF 
COEFFICIENTS MATRIX* 


BY 
E. G. KOGBETLIANTZ 
International Business Machines Corporation 


Introduction. To solve a system of linear equations Gr = g, where G is a complex 
matrix, usually an equivalent system Hx = h is formed in which H = *G’G is Hermitian 
and h = *G’g. The solution of Hx = h is performed again by replacing its complex 
equations by real ones 


Au — Bv =m, Bu + Av =n, 


where H = A + iB, x = u+ iv, hh = m + in. Thus, finally the order of G is doubled 
which is undesirable even when electronic computers are used. 

It is preferable to solve Gr = g without doubling the order of G and this can be done 
inverting the complex matrix G by diagonalization. If such two unitary matrices U, T 
can be formed, that 


*U'GT = K, 


where *U’ denotes the conjugate transpose of U and K is a diagonal matrix, then 
G”' = TK~'*U’ and the direct solution of Gz = g is obtained: 


z= Gg = TK '*U’g. 


This form is very convenient since the multiplication of matrices is performed by an 
electronic computer in almost no time. 

The unitary matrices U and T are known to exist (for example modal matrices of 
G*G@’ and *G’G), but the question how to form them for a given complex matrix G of 
large order is a difficult one. In what follows they are formed as infinite convergent 
products of simple unitary matrices u,(z,) which represent plane rotations through com- 
plex angles z, = 0, + 7, . Therefore, their practical computation is easy for electronic 
equipment and it yields a new method of solution of linear systems. 

When this method is applied to a Hermitian system Hx = h or to a symmetric 
system of real equation (particular case of Hx = h) the diagonal matrix K into which 
H is transformed has as its elements the characteristic roots (Eigenvalues) of H, while 
U = T yields the characteristic vectors, because K = *U’HU. The same holds for skew 
symmetric Hermitian matrices [(a,, + 7b,,)] with b,, = b,,,a,, + @,, = 0. 

To achieve the diagonalization of a given complex matrix A = A, a sequence [A,] 
of its transforms is constructed, the transform A; being formed from A,_, by right and 
left multiplications: 


A; *ul(¢;)A;-1u,(2;). (1) 
Thus 
A, = *U‘AT, and K = lim A, , 


with 


T =limT, = lim [[ wu; U = lim JU, = lim [] u(é,). 


n= oO n= co 1 n=o n 





*Received March 29, 1954. 
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The parameters 2; = 0; + 76; and ¢; = & + 7m; in (1) will be chosen at each step 
in such a way that K is diagonal. In practice finite products *Ux and Ty are used, N 
being sufficiently large to insure a good approximation Ay to K. This means that all 
off-diagonal elements of Ay (which tend to zero when N — ~) are already less in absolute 
value than a prescribed quantity, for instance less than 10~””. 

1. Unitary matrix U(z). Let z’ = @ — i denote the conjugate of z = 0+ 7. To 
form u,,,(z), a matrix of order n X n, we replace in the identity matrix of the same 
order ((6;;)) the four elements 6,, = 1, d:, = 0, 6,4 = 0, and 6, = 1, wherel Sk < 
m < n, by the elements 7? cos z, —y' sin z’, y' sin z and 7? cos z’ respectively of the rota- 
tion matrix through a complex angle z, the symbol y denoting the hyperbolic secant of 
2:7 ='sech 2¢. Thus, u;,,(z) is defined as follows: 





jl, | 
| 
Wee | 
| i 
ae 
9 > . 
|\— — — 7'” cosz — — — —y'”’ sin z’ — — —| kth row 
| i | 
Urm(zZ) = | i 
| ~s1 
- ——y,' 2 sin z — — — x” cos 2’ — — — mth row 
| 
i es 
| om 
~ 
| | ars | 
kth mth 
column column 
It is unitary since y(|cos z |” + | sin z |”) = 1. The matrices u;(z;) and *u}(¢;) = *u;(¢/) 


in (1) are in fact u;(z;) = Uz;m,;(2;), *uj (f; 
m, vary with j. 
Four conditions needed to determine z and ¢ are obtained as follows. Transforming 


= *u,,,,(¢/) since not only z; , ¢; but also k, , 


A = [(a,,)] into C = [(c,,)] = *uin(¢’)AUem(Z), we make vanish in C the two elements 
Cum and c,,.(k < m) and this gives four real equations with four real unknowns @, ¢, & 
and n. The sequence [Ay] is cyclic, each group of n(n — 1)/2 consecutive steps forming 


a cycle. In a cycle we make vanish in a prescribed order one after another all n(n — 1)/2 
pairs of off-diagonal elements symmetrical with respect to the principal diagonal. 

A particular pair vanishes only once during a cycle. It reappears again in the next 
step, but weakened since—as will be shown later—its modulus, when it reappears, is 
diminished in an average ratio of 1 to 2~*. At each step the sum of squares of moduli of 
diagonal elements increases, the amount of vanishing elements | a, |? + | a, |* being 
transferred into the principal diagonal. Since in our transformation the norm of A is an 
invariant, this means that at each step the sum of moduli of off-diagonal elements 
diminishes. After a certain number of cycles the moduli of all off-diagonal elements 
become less than a prescribed, arbitrarily small, but fixed, quantity. In other words, 
the off-diagonal elements tend uniformly to zero when the number N of steps increases 





without limit. 
The convergence of the sequence [Ay] to a diagonal matrix K is not sufficient to 


insure the practical importance of our diagonalization method. Should its speed depend 
on the number p of cycles already performed, decreasing when p increases, it might be 
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practically impossible to reduce the moduli of off-diagonal elements below a prescribed 
quantity. 

It is a favorable circumstance that the speed of convergence of our method does not 
depend on the number of steps already performed. On the average it remains constant 
however large the number of cycles becomes: as will be shown later, the moduli of off- 


diagonal elements decrease in the ratio of one to e~* per cycle. 
Thus after, at most, 46 cycles the average modulus of all off-diagonal elements of 
Ay , where N = 23n(n — 1), will be less than Me~** = M10~"’, M denoting the largest 


modulus of off-diagonal elements of A. 
2. Computation of C. We formulate now the rules for the numerical computation 
of elements c,, of C = *tumn(¢’)AUim(z), omitting the derivation of these rules. For any 


z and ¢ we have the following invariants: 
CikCmm — CrmCmk = AkkAmm — AkmAmk » 
l cin |? + | eim |? = | aye |? + | Gin |’; | ces |? + | ems? = | aes ? + | ns |’, (2) 
Cir ‘ + Chim , + | Conk |’ a | Cones ? _ | One |? + | Gem |’ + | Qe |’ + | Qne ¢ (3) 
Equations (2) and (3) prove the invariance of the norm: || C || = || A ||. Equation (3) 


fOr Crim = Cmk = O becomes 
2 12 | 2 12 2 2 
As | + Arm | + Amk | + { =3 = | Cex | + | Cam | 
which shows that the amount | a;,, |” + | an. |” of off-diagonal terms which vanish in C 
is transferred into the principal diagonal. Denoting the real and imaginary components 
of a complex number z by 2’ and iz’ we first form 33 real numbers in six steps. 


1) M = (a + hnm)/2, r= (ai + an,,)/2, 
P = (aj, — a,,)/2, p= (ati — av,)/2, 
N = (Qin + Gme)/2, mm = (ain + ani)/2, 
Q = Gin — Onx)/2, q = (aj, — at,,)/2. 
2) A = PM +pr—-(NQ+ nq), B* = MN +r — (PQ + 7-9), 


A* = PM + pr+ (NQ+nQ), C =rQ — qM + (pN — nP), 
B= MN +m-+ (PQ + pq), C* = rQ — qM — (pN — nP). 
3) "= (A* }. B?)' J w* _ (A*? + Bey". R - (w* + Cr _ (w*? + ofr: 


cos 26 = 6 = A/w, cos 2§ = 6* = A*/w*, 
1) 
sech 26 = y = w/R, sech 2n = y* = w*/R. 
5) sin 6 = s, = o(B)[(1 — 8)/2]'”, sin = s* = o(B*)[(1 — 6*)/2]'”, 


[(L + 8*)/2]'” > 0, 


cos 6 = c, = [(1 + 8)/2]'” > 0, cos § = c* 
(4) 


s¥ = o(C*)[(1 — y*)/2]"”, 
cf = [(1 + y*)/2]” > 0. 


o(C)[(1 — y)/2]'”, y*’’sinhn 


1/2 


(i +y/2/?>0, y*'coshn 


II 


I 


y'’sinh¢ 


y'cosh¢ = cz 
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Here the symbol o(v) denotes the sign of a real number », so that o(v) = v/| v |, if v # 0, 
while o(0) = 0 by definition. All eight numbers computed in (4) are less in absolute 
value than one. The two cosines c, , ct are positive, acute angles being chosen as @ and £. 


6) G = (M + Phect + (M — P)s,st + (N + Q)sicf + (N — Q)ast , 
g = (r+ pect + (r — p)sist + (n + Qg)sict + (n — g)c,s* , 


(5) 
H = (M — P)e,c* + (M + P)s,s* — (N — Q)sict — (N + Q)esT , 
h = (r — poec® + (r + p)sist — (n — g)sict — (n + gest . 
Now ¢x, and Cm are: 
y+ yen = 2[(G + ig)exct — (H + th)s,s3], 
(6) 
(y + y*)\emm = 2[(H + th)c.ct — (G4 + 7g)s283]. 


Among the off-diagonal elements of C only those in the kth and mth rows and columns 
are to be computed since c;; = a;; for 7,7 # k, m. For each value of 7 # k, m,1 Sj Sn, 
the elements C;; , Cjm 5 Cx; » mj are computed at a time as a group. First eight real numbers 
a, \4; a*, A*; B, uw; B*, u* are formed: 


at Ar = €,0;; + 8)A;m , a®* + ih* = cta,; + stan; , 
B+ ip = C,Ajm — 8:Ajx , B* + in* = cia,,; — sta; , 
and then 
Ci. = (a+ ade. + (8B + w)is, , C.; = (a*® + idr*)c% — (B* + ap*)ist , 


(7) 


II 


Cim = (8 + ip)e. + (a2 + A)is, , Cais (B* + ip*)c% — (a*® + 2dA*)isé . 


To form the matrices *U/ and T,, the elements of wu,,,(z), and *u,,,(¢’) should also be 


computed. They are given by 


1/2 ° ° 1/2 ° ia ° 
y’" sINZ = 8,Co + 1€)82 , _*" sin ¢ = sich + wis} , 
v'/? cos z = CiCo — 18,82 , y*'** cos ¢ = cich — i883 . 


i 
3. Exceptional cases. The rules of Sec. 2 hold in the general cases when ww* ~ 0. 
If ww* = 0, some of them are modified as indicated below. We discuss five exceptional 
cases: 
I) w* = 0, w + 0; ID) w = 0 u* € 0; IID) w = w* = 0,C = C* = 0;IV) w = w* = 
C = C* =0, butn’ +Q@>0;V)w= wt =C=C¥=n=Q=0. 


Case I. | C* | = R > w > O, so that o(C*) has a meaning. We take s, , ¢, , 82 , C, as 
in (4), but s§ = o(C*), 2) ct = 1/2’ since y* = 0. Now B* = 0 and o(B*) is meaningless, 
but in this case £ = 0, so that s{ = 0 and ct = 1. The elements c;,, , c;,, are given by 
(7), but 

‘ 1/2 ’ ‘ 2 Y 
ch; = 2°" lai; + atio(C*)], eff = 2°-'*[ali — aljo(C*)], 


(8) 


Cny = 2°" [an; + atfo(C*)], ems = 2°" [an — atjo(C*)], 
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as well as 


veh, =.2'”[Ge. — Hs.0(C*)], yet, = 2"*[9c. — hs,o(C*)], 


ll 


yCm = 2’? [He. — Gs,o(C*)], yell, = 2'[he, — gs,0(C*)]. 
Case II. This is similar to J:y = Oand|C| = R > w* > 0. 
We take s*, ct , s€ , cf asin (4), but s, = 0,c, = 1, s = o(C)/2' and c, = 1/2*. The 
expressions of elements C;, , Cim » Cix » Cmm ATC: 
ci, = 2°" la}, — aj/,o(C)], ef = 2°" [ali + a}na(C)], 
f= 2°" Ia! — alic(C)], cf, = 2°" fall, + a},o(C)], 
y*ch, = 2'”*[Gct — Hsto(C)], v*Ciun = 2’ [gC. — hste(C)], 
7*C,m = 2’ [Hct — Gsto(C)], Y*Cnm = 2' [hes — gsto(C)]. 
Case III. Here y = y* = Oandc, = ct = 1, 8, = st = 0, ce. = ct = 1/2!, while 
8. a(C)/2', st = o(C*)/2*. The expressions of diagonal elements ¢,, , Cnm Te: 
ec, = M+P4+ (M — P)o(CC*) — (n+ go(C) + (n — go(C*), 
ck =r+pt (r — pio(CC*) + (N + Qo(C) — (N — Qo(C*), 
ce, = M—P+(M + P)ao(CC*) — (n — ™ol(C) + (rn + Do(C*), 
cf =r — pt (r+ pio(CC*) + (N — Qa(C) — (N + OAo(C*). 
The off-diagonal elements are computed as in (8) and (9). 

Case IV. The conditions c,, = 0, Cn, = 0 yield one and the same equation, so that we 
add the condition ¢ = z. We subdivide this case n? + Q’ > 0 into: IVa) n ¥ 0 and 
IVb) n = 0, but Q ¥ 0. 

In the first case [Va): 
t= b= pet py; raya +p), 
where R = (n® + p? + Q’)*. Formulae (4) hold for c; , c2 , | s; | and | s2 |, the signs of 
s, and s, being now o(n) and o(Q) respectively. If Q = 0, then s, = 0. The off-diagonal 
terms in C are computed as in (7), while 


af M + NR/n: ci=r+R; ce, = M — NR/n; cl. =r—R. (10) 


C.. = 
IVb). Now n = 0, Q ¥ 0. Interpreting o(0) as zero, we have 6 = o(p) for any p 2 0, 
so that 


st = 8, = {[1 — o(p)]/2}'”, ct =e, = {[1 + o(p)]/2}™, 


while cf = c. = [(1 + | p |/R)/2]}, st = s2 = o(Q)[(1 — | p |/R)/2]*. The elements of 
C are computed as in IVa), replacing in (10) N/n by — ¢/Q, that is: 


7 


ch =M—QR/Q; ch=r+R; chan =M+QR/Q; C=r—R. 


Case V. In this case the vanishing of n and Q entails also p = M = Osincew = w* = 
C = C* = 0. Thus, c,, = 0 and c,, = 0 reduce to only one equation, namely 


P sin 2z — N cos 2z + iq = 0, (11) 
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if we add the condition ¢ = z. This condition ¢ = z presupposes N* + P” > 0, since if 
N = P = 0, then g cannot vanish because of | an, |? + | dem |? = N? + 4° > 0. 

Va). If w* = N’? + P’ > 0, then, denoting w* + q’ by R’, we have: 6 = P/w, y = w/R, 


a(sin 8) = o(N), o(sinhy) = — o(q). Since 6 = £, g = 7 we have ct = ¢, etc., andc, = 
[1 + 6)/2]*; s, = o(N)[1 — 8)/2]*; cg = [(1 + y)/2]*; 82 = — o(g)[1 — y)/2]'. Now 
G= —H=w,g = h =rand thus cy, = R + 7r, whilec,, = — R + ir. 


Vb). If N = P = 0, we take ¢ + 2’ = O that is¢ = n, 0 + — = O. The condition 
Com = Cme = 0 becomes tan 26 = qg/r so that g = 7» is arbitrary. We take g = n = 0. 
Since 6 = r(r? + q’)~? we have c¥ = c, = [(1 + 6)/2]}, st = 8, = o(q)[1 — 8)/2]!. The 
elements of C are: Ci: = Cam = O(r? + Q°); cj, = @ + 1A Cim = BH tu; = a* + 2A* 
and c,,; = B* + i*. 

4. Convergence. Let us consider two consecutive transformations of A into C with 
C.; = C;, = 0 and then of C into D with d,,, = d,, = 0. Thus, the elements a,,, , Qnx , 
Qi; , Um; » 4; » Ajm Of A are transformed according to the scheme: 


a; ~0— d,; ; Aim — Crm > 0; Aim > Cim 2 dim 3 
a;,70—- di ; Omk —> Cm~e — 0; Qmj —> Cmj —> Am; - 
The relations (2) yield | dj, |" + | dim |" = ¢im |"; | de; |) + |] da; | = | en; [°. Now 


2 oa 
c = 7* | a;,, cos ¢’ + a,,, sin ¢’ |", 


im 


Cm; | = ¥ | Qn; COSZ + a,, SiN z 
Using polar forms so that a,, = | a,, | e'°"* etc., we have: 
2 | (2 » 74 | P49 *] cos 26 
| Cmij = ani mk Qmj | Omk COS 2 
(12) 
+ 2 | GniQne | [y sin 26 cos (&n; — War) + tanh2yg sin (wn; — Wms) ] 


and a similar expression for 2 | ¢;,, |” in which a,,; , a,,, , 9 and g are replaced by @;n, Qim , 
£and — 7 respectively. Each pair of symmetric elements a,, , a,, , $ # 7, vanishes and 
reappears once in a cycle. When the number of cycles p increases the various values of 
6, v, £, 7, re , We, related to a fixed pair a,, , a,, can be considered as random values. 
Therefore, when p increases the average values of sine and cosine of 26, 2, w,, — w,; , 
w,, — w;, and of tanh 2¢ tend to zero. It means that when p is large the two last terms in 
(12) and in the corresponding expression 2 | c;,, |" do not contribute to final values of 
moduli and these values are to be derived from the approximate expressions 


| crm |? = [| @im |? + | Gem |7)/2; | Gms |? = [| ans |? + | Gua |7J/2 
| ] 


which represent the average effect of one cycle, when a large number of cycles is con- 
sidered. This gives, for one cycle, 


| dix |? + | dim F + | d,; |? + | Oni ’ - [ Aim . so | Aim ; + | ami |? + | Ami °] 2. 


Consider now the global effect of a cycle on the average modulus » of all n(n — 1) off- 
diagonal elements, where 


v? = (>> | a,; |)/n(m — 1). 


+) 


Since for each pair of symmetric elements the sum of four moduli squared of elements 
of D is approximately equal to half the sum of four moduli squared of elements of A, 
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the average modulus » is multiplied (1/2)n(n — 1) times per cycle by the factor [1 — 
2/(n — 1)nj'. 

At each step, indeed, two among n(n — 1) squares of moduli of off-diagonal terms are 
lost, the sum of four new ones being equal to the sum of only two old ones on the average. 
After the completion of a cycle with its n(n — 1)/2 steps, the average modulus v becomes 


v[l — 2/(n — Inf?“ wove”, 


Thus, after p cycles the average modulus of off-diagonal terms decreases e”” times. The 
speed of convergence as characterized by a constant factor of decrease, namely e* per 
cycle, does not depend on the number of cycles. The results of numerical computations 
performed with the aid of IBMs new electronic data processing machine type 701 agree 
completely with the conclusions of this study of convergence. The type 701 performs 
with extreme rapidity: a 32 X 32 matrix is diagonalized in 19 minutes and the numerical 
results are printed in four minutes, the total time being 23 minutes. Since the errors do 
not accumulate the procedure is very accurate: reconstructing the original matrix from 
the diagonal matrix and the Eigenvectors, if the matrix is normal, or from the diagonal 
matrix and unitary matrices U and T in the general case, very small differences are 
obtained between the elements of the original and the corresponding elements of the 
reconstructed matrix. Thus, for instance, in the case of a real symmetric matrix of order 
32 x 32 these differences were less than 10~°. Since the elements of the original matrix 
were of the order of 10~° the relative differences were less than 107’. 

5. Special matrices. Our diagonalization method does not preserve the skew symme- 
try. A skew symmetric matrix is diagonalized as a general complex matrix, its skew 
symmetry being lost from the first step. Since for a skew symmetric matrix a: + Qnx = 0, 
a,, = 0 its transformation is characterized by 


Ch; + Ci, = Ug; + 04,,; 5 Cui TCim = —v'; + Udy; , where 


u = y* cos ¢’ — y' coszand v = y*! sin ¢’ — y' sin z. If we want the skew symmetry 
preserved, we must have ua,; + va,,; = 0 and v’a,; — u’a,; = 0 for any a; , a,,; So that 
u = v = 0 and this gives y* = y, ¢’ = z. But then ¢’ = z entails together with a, = 
Qnm = 0 the corollary — ¢,,, = Csm = @m and we cannot make vanish in C the elements 
Cmte» Cem » But the symmetry of a complex matrix as well as the Hermitian and skew 
Hermitian character of matrices are preserved. We consider in this paragraph the 
symmetric complex, general real, skew Hermitian and Hermitian matrices. In each case 
the general procedure is somewhat simplified thanks to special properties of the class 
considered. 

Symmetric complex matrices. If dim = Anz , then Q = g = 0,C = C* = 0, w = w*. 
Therefore ¢ = 2’ and ct = c, , st = 8, ,c$ = C2, 8% = — 8, which entails also a* = a, 


B* = B, \* = \ and yw* = u. Formulae (7) hold and they prove that c;, = c,; . Formulae 
(6) are simplified: 


Cur = Chr + UK, fn he + > with 
cl, = P+ (RM + DN)/w, cl, = —8P + (RM — DN)/w, 
cl, = 6p + (Rr + Dn)/u, ch’, = —dp + (Rr — Dn)/w. 


For real general matrices the numbers n, p, g, r vanish, so that C = C* = » = » = 0, 


6¢=& co = ct = 1,8, = s§ = 0, wo = w* = R = R,R, , where 


y=7* =1,2 
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Ri = M’ + Q and R? = N’ + P’. We have also 6 = cos 20 = (MP — NQ)/R,R, ; 
6* = cos 2 = (MP + NQ)/R,R, as well as 
8, = o(s,)[(1 — 6)/2]'”; ¢, = o(,)[(1 + 6)/2]'”; 
(13) 
st = o(s*)[(1 — 6*) i c* aa a(c*)((1 + 5*) ay”. 


? 


The signs in (13) must agree with the relations M sin (@ — £) = Q cos (@ — £) = MQ/R, ; 
P sin (6 + £) = N cos (0+ £) = NP/R, and thus, for | MN | # | PQ |, the signs are: 





a(s;) a(s*) a(c;) a(c*) 

For MN |>|PQ|>0: (MN) o(N) | +1 | o(M) ae 
For PQ|>|MN|>0: o(PQ) —o(Q) +1 | o(P) 
Bos Ma, Maia 

Suppose | MN | = | PQ | > 0, the singular case MN = PQ = 0 being discussed later. 

If | MN | = | PQ | > O then both rules (14) hold since this is a limiting case when 


| MN |—| PQ|for| MN | $ | PQ |. No contradiction can arise between the two rules 
(14): o(st) becomes meaningless, if MN = PQ, because then sin = s*¥ vanishes, while 
for MN .= — PQ, s, = sin 6 vanishes so that o(s,) is meaningless. 

In the singular case MN = PQ = 0 the parameters N and Q cannot vanish simul- 
taneously because 2(N* + Q’) = aim + an, > 0, otherwise no transformation is applied 
tO dim = @,, = O and this pair is skipped. At least two of four parameters M, N, P, Q 
vanish since MN = PQ = 0, but it may also happen that three of them vanish in which 
case either VN ~ 0, or Q = 0. In all we distinguish five subcases: 


1) Ait = Onm » Gem + Ont = 0, 80 that N = P = 0, but MQ = 0; 
2) Giz + Anm = 0, Gem = Anz , 80 that M = Q = 0, but NP = 0; 
3) Giz = AOnm = 0, Gem + Anz ~ 0, So that M = P = 0, but NQ ¥ 0; 
4) Que = Onm = 0, Qim + Gunz = 0,80 that M = P= N =0,Q #0; 
5) Qiu = Gan = 0, Gi, — 222 = 0,80 that M = P=Q=0,N #0. 


In all these five cases st = 0 and c¥ = 1, while the values of s, = sin @ and c, = cos 6 
are found to be: 








Case: 1) 2) 3) 4) | 5) | 
8, = Q/R, | N/R, o(Q) o(Q) | o(N) | 
—_ —EE — = a 7 —EEE ————E —— 7 _ 

} 
¢, = M/R, P/R; 0 | 0 | 0 


The elements of C are computed using the following expressions: 
Cie = C105 + Sidim ; Cin = —80;; + CiAjm. 5 


Cr; = Cid; + 81m; , Cui = SIO; + CIO; , 


(M? +. Q’)' 2 + (N? + 20 a _ (M? + ey" ws (N? + P*)' a 


° 
ks 
ll 
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If a real matrix is symmetric, then the procedure is simplified further and for a 
real symmetric matrix we have 6 = P(N? + P’)~* so that 

ct = c¢, = [(1 + 5)/2]'”; st = 8, = o(N)[(1 — 8)/2]'”. 
The elements of C are as follows: 


M + (N’ + PP”), Cum = M — (N’? + P”)™”, 


II 


Cir = C10jx + 8:Ajm , Cim = —8:0;, + CiAjm « 

The diagonalization yields the characteristic roots and vectors of the original matrix 
since u’(£) = u’(@) = [u(@)]"’. 

Skew Hermitian matrices. The skew symmetry of a Hermitian matrix is preserved. 
In this case M = P = N = q = Oand ay = i(r + Pp), Qnm = U(r — P), im = QD + ™M 
and a,,, = — Q + in, with Q’ + n’ > 0. Here we find H = G = 0 so that c,, and Cam 
are purely imaginary. We take ¢ = z since the vanishing of ¢,,, entails also c,,, = 0. It is 
easy to verify that a* + a = 0, 8* + B = 0,A* = A, uw* = wso that, +c, = 0 and 
0. This completes the proof that the skew symmetry is preserved. Since 
¢ = z we have *u(¢’) = u ‘(z) so that our procedure yields the characteristic roots and 
vectors of skew symmetric Hermitian matrices. It proceeds as usual in the case of a 
complex matrix, a half only of off-diagonal elements being computed because of the 


e T Cin 


™ 


skew symmetry, while the diagonal elements are: c,, = i(r + R) and Cam = i(r — R) 
a “Seu 2 2 2\4 
with R = (n° +p + QQ)’. 

Hermitian matrics. Now a,,, = Qim aNd 4,4 , Gmm are real. Therefore, m = n = p = 

Q = 0 and ¢ = z since the vanishing of ¢,,, entails that of c,,, = Cim , the character of 

a Hermitian matrix being preserved. We have, indeed, in this case a* = a, B* = 8B, 


A* +r =0,u* + yu = Oso that c,; = ci, andc,; = cim.On the other hand, g = h = 0 
and Cx: , Cmm are real. Denoting (N? + P*)' and (N? + P? + q°)' by w and R respectively, 
it is found that 6 = P/wand y = w/R which gives ct = c, , s§ = 8, , c$ = Cc. and s} = 8, 
in their usual form with o(s,) = o(N) and o(s.) = — o(q), while c, > 0, c,. > 0. These 
results hold in the general case, when N # 0. When q = 0, we have s. = 0 so that o(q) 
is not needed. The elements of C are computed by (7), while the diagonal elements are 
simply c,. = M + R, Cnn = M — R. 

In the case when NV = 0, we have q ¥ 0 because | ai» * = | On: |? =N’+ q cannot 
vanish. The numbers ct = c, and st = s, are given by the formulas s, = — o(q)[(1 — 
| P |/R)2]', ce. = [1 + | P |/R)/2]}' for P 2 0, while st = s, = [1 — o(P)]/2 and ct = 
ce, = {1 + o(P)]/2 with c(0) = 0, as usual. The diagonal elements are again c,, = M + R, 
= - M — R where R = (p’ + q’)' since N = 0. The off-diagonal elements are com- 
puted as in (7) in which this time, if P # 0 


Qa + ir) = ay + Asm + (a;, — A;n)0(P), 
(P ~ 0) 
2(8 + ip) oe + Aim + (Qin + djm)o(P). 


For P = 0 we have (a + id)2' = a;, + a;, and (8 + ip)2* = ain — Qj. 
6. It is interesting to observe that the diagonalization of the most general complex 
matrix @ yields the upper and lower bounds for the absolute values | \; | of its characteris- 





tic roots. 
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A known theorem (H. Weyl, Proc. Natl. Acad. Sci. U.S. 1949, pp. 408-411) states 


that the common characteristic roots a; of two matrices G*’G and GG*’, aj = a; = 
a; > -+: = a& = O are related to those d, , \., --- A, of G, [A.J] 2/aA.|] 2 /As| 2 
> |A,|, by an equality []7|a,; | = [J]? a; and by n — 1 inequalities ir ikhis 

2) 


—s 


Il? a; (m =n — 1,n — 2, --- 3, 2, 1), where a, denotes the positive square root (a; 


Therefore, a, and a, are the bounds for the moduli of X,; : 


IV 


a, 2 Mi/2°::2 An | An - 


On the other hand it is easy to prove that the moduli of elements k; (these elements 
are complex numbers in general) of the diagonal matrix K into which our diagonalization 
transforms G are precisely equal to a; :| k; | = a; . 

Thus, although the diagonalization in general does not yield the characteristic roots 
A; of G, it gives the upper and lower bounds for their moduli: | k,| S|, | S|, | S | &, |. 

All the proofs were omitted for the sake of brevity and also because of their elementary 


character. 
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INDENTATION PRESSURE OF A SMOOTH CIRCULAR PUNCH* 


BY 
E. LEVIN 
National Bureau of Standards, Los Angeles, Calif. 


Abstract. A smooth, flat, rigid punch under increasing normal load presses against 
a half space of a perfectly plastic material which obeys Tresca’s yield criterion. An 
admissible velocity field is constructed for an arbitrary smooth punch hence, for any 
particular case, a limit design theorem of Drucker, Prager and Greenberg may be used 
to compute an upper bound for the punch indentation pressure. A lower bound for any 
convex area of indentation has been given by Shield and Drucker. The results of the 
present paper are used to compute an upper bound for a punch with circular cross 
section. It is conjectured that this is an upper bound for an arbitrary punch. 

1. Introduction. The exact solutions of any but the most trivial problems involving 
plastic deformation have thus far been impossible to obtain. Consequently the limit 
design theorems of Drucker, Prager and Greenberg [1] which provide upper and lower 
bounds on the collapse load have become the most powerful tools available for the 
analysis of such problems. In the present paper this technique is used to obtain an upper 
bound on the punch pressure at the moment of impending plastic indentation. 

Consider a flat rigid punch subjected to increasing normal load and pressing against 
a half-space (Fig. 1). The half-space is assumed to be of an elastic-perfectly plastic 











Fia. 1. Punch on half space. 


material which obeys Tresca’s [2] yield criterion of constant maximum shearing stress, 
k, during plastic deformation. The pressure is gradually increased until indentation due 
to plastic flow of the material is impending. Geometry changes up to this point are 


*Received April 10, 1954. The preparation of this paper was sponsored in part by the Office of Naval 
Research, USN. 
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presumably small and are neglected. The actual load applied to the end of the punch 
may be concentrated or distributed but it is assumed that over the contact area there is a 
uniformly distributed normal pressure. 

According to the Drucker-Prager-Greenberg theorem an upper bound on this inden- 
tation pressure, p, may be obtained by constructing a kinematically admissible velocity 
field. Such a field must satisfy the following conditions: 

(a) it must accommodate the assumed downward motion of the punch; 

(b) it must be incompressible in the sense that 


6 +téea+te = Q, (1) 
where ¢;, €, €; are three orthogonal strain rates; 
(ce) | pv dS > | 2k max |e | dV + | k Av dSp , (2) 
’§ Pe | ’“ Sp 


where S is the contact surface between the punch and the material, V is the volume of 
the deforming material, Sp represents surfaces across which the tangential velocity 
is discontinuous, Av is the magnitude of the relative change in velocity across Sp, max | ¢ | 
represents the absolutely largest principal strain rate component. 

Actually the third condition is used to compute p when a velocity field satisfying 
the first two conditions has been found. If this velocity field involves slip between the 
punch and the material, it is necessary to make the further restriction that the punch 
be smooth. 

2. General case. In this section an admissible velocity field is constructed for an 
arbitrary punch. Since the velocity field involves slip between the punch and the material 
it is assumed that the punch is smooth. The only restriction on the shape of the punch 
is that its cross-section be star-shaped (Fig. 2). Thus the boundary of the cross-section 


(a) (b) 


Fic. 2. (a) Admissible cross section 


(b) In idmissible cross section. 


can be described in polar coordinates by r = d(@) where d is a continuous single-valued 
function. This is a very weak restriction and admits, for example, any punch of convex 
cross section. 

For a given punch a suitable origin is selected and a cylindrical coordinate system is 
established at this point. In each section @ an assumption is made as to the general 
shape of the velocity field (Fig. 3). The exact form is then determined so as to satisfy 
the requirements for a kinematically admissible velocity field. 


In region A assume 
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Then to satisfy incompressibility 


OF c oF _ 


ort 0z 


or 
A solution is 


F = ¢(r — 2)/r. 
To accommodate the assumed downward motion of the punch [V,],.o = v. Hence 
F = v(r — 2)/r. 


In each of the remaining two regions a similar assumption as to the general shape 











~ 


Fic. 3. Velocity Field. 


of the velocity field is made in accord with Fig. 3. The exact form is then determined so 
as to satisfy the incompressibility condition. The results are presented in Table 1. 





TABLE 1 
Region V. V, | We 
A F F | 0 
| 
B 2G (d — r)G 0 
9! _* om 9! 25 we 2 y 1/2 
where F = v(r —2)/r, G=- me — 3 ie 7 Sh } H = v(2d — r —2)/r. 


2 2)1/2 
r{(d — r) z} 


As well as satisfying the incompressibility condition, it can be verified that this field 
is continuous between the three regions and vanishes on the boundaries r = z, 





136 E. LEVIN [Vol. XIII, No. 2 
td — r)? + 2]? = d/2'” and r + z = 2d. Hence the last integral in Inequality 2 
vanishes and an upper bound on the indentation pressure is furnished by any pressure, 
p, satisfying 


ee | [d(@)}? do > 2k | max |e | dV. (3) 


3. Circular punch. For the particular case of a punch with circular cross section, 
P(8@) is a constant which will be taken as unity. The strain rates appear in Table 2. 


TABLE 2 


= ar — = : a ee an ites 
{ B ( | 
-_ | 
| 1/2 2 ‘ } | 
2 2Qvz 2° v2(2r — 8r +142) v(2 — 2) | 
€,, —-— 7 a 
| °° rid —r? +2} °° 
| 
OT — 2 a. on 21 AZ— Tr — Zz) 
| €9¢9 — > 2 3)1/2 ~ 3 2 
r rill —r) +2} r r 
| 7 2 vz( — 7) ? 
« - — ‘ 
ri(l—r) +2z);" Y 
| 2 1) 1/2 »! ) 
ge=—fr) wl —?)., 2°v1 — 7 ils W2—r-—2z 
€, >» Was fase 2 1 4238/2 —iC@ ae Oe : yp" 
Ty \ T) — ot , oho SS in a a! 


The principal strain rates are 


€ —€ 
é& = [—€ps + (e,, + 4e, — 4e,,€ )’ |, 
af < i 2 4 1/2 
e, = 4[—€99 — (eso + 4e,, — 4e,,€..) J. 
Since ¢,, > 0 in each region, the absolutely largest principal strain rate will be «; pro- 
vided 4¢°, — 4e,,€,, > 0. This is easily verified in regions A and C. In region B it is found 
that 
2 0G 0G 0G 0G 
fe, — 4e,,€6,, = | 2— + —7)5 — 421 — r) — = 
Oz Or Or oz 
0G 0G |° 
=f2—--— (1-7 >< > Dp. 
Oz or ~ 
Hence in all three regions 
Max |e | = | e, | = 4feon + (eso + 4e-, — ae iy 8 


An upper bound is obtained as any value p satisfying 


pv a 
c 2r > 2k-2r | e, |r dr dz. 


aad 7A,B,C 
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A portion of this integration was performed numerically. It was found that 5.84k is an 
upper bound on the indentation pressure. 

4. Discussion and conclusions. The assumed shape of the velocity field is patterned 
after the field constructed by Hill [3] for the two dimensional punch. Due to symmetry 
this field is most easily applied to a circular punch. For other cross-sections the general 
formulation given in Sec. 2 would involve considerable labour and indeed might not 
yield the sharpest upper bound. 

For the two dimensional punch it is known that the indentation pressure is (2 + )k. 
In a recent paper by Shield and Drucker [4] on punches with rectangular cross section, 
a lower bound of 5k was obtained for any convex area of indentation. Upper bounds 
varied between (2 + )k for a very long rectangle (which agrees with the two dimensional 
solution) and 5.71k for a square. 

It may be noted that the results of Shield and Drucker indicate that the upper bound 
increases as the ratio of cross-sectional area to perimeter increases. Since the circular 
cross-section has the greatest area for a given perimeter, the upper bound of 5.84k for 
the circle appears quite reasonable as compared with the upper bound of 5.71k for the 
square. Furthermore, it is conjectured that this value for the circle is an upper bound 
for an arbitrary punch. Combining the results of Shield and Drucker with the results 
and conjecture of the present paper indicates that the indentation pressure for any 
smooth punch with convex cross-section lies between 5k and 5.84k. 
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BOOK REVIEWS 


Probability and information theory with applications to radar. By P. M. Woodward. 
McGraw-Hill Book Co., Inc., New York, and Pergamon Press, London, 1953. 
x + 128 pp. $4.50. 


This monograph is generally descriptive and demands no special mathematical preparation for any- 
one acquainted with the field. The first part reviews lightly the elements of probability and waveform 
analysis that have application in radar. There follows a chapter on information theory; however, for a 
deeper understanding of this theory the reader may well want to refer to the original papers by Shannon 
and Weaver which are cited. Succeeding chapters discuss: 1) the statistical problem of reception, 2) a 
simple theory for range measurement of a stationary point target in the presence of white Gaussian noise, 
3) the mathematical analysis of radar information and 4) some aspects of the transmitted radar signal. 

The book satisfies its goal of illustrating the way in which probability applies to communication and 
radar theory without requiring a highly advanced mathematical background but should not be regarded 
as a fundamental text for a student in the field of communication and radar. 


L. C. MaxImMon 


The mechanism of economic systems. By Arnold Tustin. Harvard University Press, 
Cambridge, Mass., 1953. xi + 161 pp. $5.00. 


Professor Tustin’s book is well described by its subtitle: “An Approach to the Problem of Economic 
Stabilisation From the Point of View of Control-System Engineering.” In his opening chapters the author 
translates into the language of systems diagrams a number of the most important theories advanced by 
modern economists for the explanation of economic fluctuations. These chapters serve the dual function 
of acquainting the economist with the techniques used by engineers to diagram dynamic systems, and 
acquainting the engineer with the general dynamic character of an economic system. The engineer who 
masters them will have advanced a long way toward an understanding of modern theories of the business 
cycle. The economist, to achieve an equivalent sophistication in modern control-system theory, will 
elementary, original in 





need to work carefully through Chapter 3; which provides an introduction 
presentation, but compact—of methods of systems analysis. 

With these chapters as foundation, the last half of the book carries both economist and engineer 
through a consideration of possible elaborations of the economic model to achieve adequate realism, and 
of the ways in which analogue computers might be used for the study of the inevitably complex and 
non-linear systems that would emerge from this elaboration. Professor Tustin makes a strong case for 
analogue computers as useful, and perhaps indispensable, tools for extending our knowledge of the 
economic mechanism to the point where we can evolve appropriate policy measures for its stabilisation. 

Although Professor Tustin is very modest about his knowledge of economics, he shows considerable 
sophistication in that field, and a wide acquaintance with the relevant literature. Hence, his book can be 
recommended strongly to the engineer and natural scientist as a route, starting from thoroughly familiar 
territory, that will carry him into economic theory; and to the economist as a starting point in acquiring 
a knowledge of the techniques of servomechanism analysis and of computer possibilities. Both of them 
will be aided in their future study by an excellent bibliography. Readers, whatever their background, 
will also find in Professor Tustin’s pages a rich source of original proposals about business cycle theory, 
which—whether in individual cases they turn out to be right or wrong—will provide many ideas for 
further investigation. 


HERBERT A. SIMON 


(Continued on p. 160) 
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THE DESIGN OF TWO-DIMENSIONAL AEROFOILS WITH MIXED 
BOUNDARY CONDITIONS* 


BY 
L. C. WOODS 
Sydney University, Australia 


Summary. A method is given of designing two-dimensional symmetrical aerofoils in 
subsonic compressible flow when the boundary conditions are mixed. The aerofoil 
surface is divided into three parts such that the shape is prescribed in the first and third 
and the pressure distribution is specified in the second. The method has the advantage 
over the usual method of aerofoil design (pressure distribution given over the whole 
aerofoil surface) of giving the designer direct control of the nose shape and the trailing 
edge angle. 

1. Introduction. Lighthill’s paper [1] on aerofoil design in incompressible flow and a 
subsequent extension of the method by the author [2] to compressible flow deal with the 
case when the pressure distribution is specified over the whole of the aerofoil. Now 
while the designer wants to control his pressure distribution over the front half or two- 
thirds of the chord he is usually more concerned about the shape of the aerofoil towards 
the trailing edge than the corresponding theoretical pressure distribution, which in 
any case is not attained in practice because of the action of viscosity. Indeed a number 
of aerofoils which have been designed by Lighthill’s method have been subsequently 
improved by modifying the shape near the trailing edge. It thus appears that a theory 
which permits the specification of the pressure distribution over the front part of the 
chord and the shape over the rear part would be of some practical value. The theory 
given below is a little more general than this as it also allows the shape of the aerofoil 
nose to be specified. The chord is thus divided into three parts in which shape, pressure 
and shape are specified respectively. The design problem is the determination of the 
shape of the middle section to give the specified pressure distribution. 

The theory is developed for a symmetrical aerofoil at zero incidence moving at sub- 
sonic speeds. It is exact for both an incompressible gas and a ‘Karman-Tsien tangent 
gas’, and hence for a real gas it is an approximate theory, applicable to most of the 
subsonic Mach number range. The basic mathematical theory is as follows. 


The equipotentials (6 = constant) and streamlines (¥ = constant) of an inviscid 
compressible flow are defined by 
do = qdn, dy = £ q dn, (1) 
Po 


where p, po are the local and stagnation densities respectively; n, s are distances measured 
normal to and along a streamline respectively; and q is the velocity magnitude. Values 
in the undisturbed stream at infinity will be denoted by a subscript ‘’. The independent 
variables of the theory are @, the flow direction measured from the direction of flow at 
infinity (i.e. 6. = 0), and Q, defined by 


a=-| "dq, (2) 


where U = g.., 8 = (1 — M”*)"”, and M is the local Mach number. 





*Received May 17, 1954. This work was performed while the author was a member of the New 
Zealand Scientific Defense Corps, attached to the Aerodynamics Division of the N.P.L., England. 
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Consider an imaginary gas for which the number m defined by 


m = Bpo/p, (3) 


is constant, and equal to the value it would have for an ideal gas in the undisturbed 
stream [6]. Such a gas can be shown to have two important properties. Firstly the com- 
plex function, 
f=2+ 7280, (4) 
is easily shown [8] to be an analytic function of 
w=o+imey, (5) 


where m. denotes the value of m at infinity in an deal gas. Secondly, if p is the gas 
pressure, the (p, 1/p) curve is a straight line tangential to the curve for an ideal gas 
at (p., 1/p.). This second property is quite well-known, but it is convenient to outline 
the proof briefly here. 

From Eq. (8) it follows that the imaginary gas is defined by 


ap Pol ’ 


where a is the local speed of sound. The tangent to the adiabatic curve of an ideal gas 
is dp/dp = a’, so replacing a’ in Eq. (6) by this gradient ensures that the adiabatic 
curves of the imaginary and ideal gases can touch. From Eq. (6) and Bernoulli’s equa- 


tion, namely 


| 
P41 gdq=0, 
p 
it then follows that 
o@ fg 
pF = AaPo , 
dp 
or 
pa = asp. . (7) 
Further integration yields 
» (2 L) 
D— Do = —G2fo\- ~— —], 
Z Z > a 


which is a straight line in the (p, 1/p) diagram, tangential at (p., 1/p.) to the ideal 


gas curve. 
The relation between Q and q is required in the following theory. From (6) and (7) 


it follows that 


a’ — q° = ac al qa ? 
and hence 
a8 = a8. 
Thus from g = aM, 
M B. (1 — p)'” B. 
tne aay sl (8) 


U BM 8 M..’ 
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(It should be noticed from Eq. (8) that a tangent gas can never have supersonic patches 
embedded in it.) From Eqs. (2) and (8) it follows after some algebra that 


‘ = sinh e cosech (Q + 6), (9) 


= sinh’ '(8./M.). (10) 


on 


If the flow is incompressible, 2 = log (U/q), m. = 1, and 


=m (2) =m) 
f= tog (8 e'’) = log( U ho)” (11) 


where z = x + iy is the physical plane. In incompressible flow z = z(w), and so (11) 


yields immediately that f = f(w). 

Since f is an analytic function it satisfies 

V's = 0, (12) 

in the w-plane, or in any plane derived from the w-plane by a conformal transformation. 
The solution of this equation with the appropriate boundary conditions follows in the 
next section. 

2. Mathematical theory. Figure 1 shows half of a symmetrical aerofoil at zero inci- 
dence. The shape of the aerofoil is assumed known in the intervals BC and DE, while 




















z-plane 
er amneamboumoampe be, 
Aco B E Foo 
w- plane 
wa 
-@ a | a b 
Aco BC D E Foo 
(6+i&)-plane 
itbiallinsiidiiiaiasete inal tele S27 ___ 
& FA © 
FIG.1 


the velocity distribution is known, or can be deduced, as a function of s, the surface 
distance along CD. The equation 


$(s) -/ q ds, 


enables g(@), and hence from Eq. (9), 2(¢), to be deduced in —a < ¢ < a, but 6) 
cannot be determined exactly in —e < ¢ < —a,a < ¢ < b, until q is known in these 
intervals. Ignoring this difficulty for the moment, we shall proceed to find the solution 
of Eq. (12) in the w-plane, with the mixed boundary conditions: 6(¢) known in —» < 
@< —-a,a<@< ©;Q(¢) knownin —-a<¢ <a. 
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The transformation 


w-— a 


: 1 wta 
5+ it = —= log |———_]}, (13) 
where 6 and £é are real, transforms the upper half of the w-plane into the infinite strip 
—o <§< 0,0 < &— < 4a. The (6 + 72¢)-plane is shown in Fig. 1, and it will be noticed 
that in this plane the boundary conditions are particularly simple. The solution of (12) 
in this case has been obtained in an earlier paper [3]; it is 


(6, == |  {0(5*) cosech (6* — 5 — it) + 2(5*) sech (5* — 6 — ié)} dé*, (14) 


where 6(5*) is the value of 6 on — = O, and 2(6*) is the value of 2 on = 47. Integrating 
(14) by parts, we find the alternative form 


> re 
f(6,6 = 2, + 16. + = | tanh”’ exp (6* — 6 — it) d6(*) 
T J 5 x 


9) ps 
—=| tan™' exp (6* — 6 — it) dQ(8*), (15) 
@ §*#= -) 
where 2. and @, are the values of 2 and 6 at 6 = ©,i.e. at@ = —a. From Eas. (13) 


and (15) 


os ((w - a)(o* — a)\' , 
te | ) tanh" 5 re. samoall f dd(o* 


3) . a 
f (a) = Q a 10 4 = (| 
ge a\J, \(w — a)(o* + a) 


“ o*=—e 
na ( \ %\) 1/2 
2 , (w+ alla — o*)\' 

a . j 


——a——oF 86S"), (16) 
\(w — aja + 9"*)) 


H* 
Q 


where, since 6(¢*) = 0 outside (b, —e), the range of integration has been reduced. Equa- 
tion (16) is the required solution of (12) in the w-plane. 

Two important auxiliary equations follow from (16) by expanding in a power series 
in (1/w). It is found that 


9 a— a \1/2 
P j 2 ’ Oo Tra 
) = 92, 4 tanl z ——— ) 16(p* 
f(w) 224 2 | inh eae d6(o*) 


cf "i esa ") 4 7 ae ga ; (# = “y" 
; é ‘ r 12 * € und *\ > 
+ l.. | tan (: + o* aX") + ir | | tanh o* +a dO( ) 


alt § Je 
iti ‘ 1/2 
] .( [\e = 2) 
—— , eas ( e ( ‘ 
-j- 7 ss + J, /\o* +a o* + a) dA(¢g*) 
anil / oe - | > 
é (¢ =) (@* + a) dQG*)? + O(1/w’). an 
Jgee-a \a + *, ) 


Since @.. = 0, g. = U, it follows from (2) that f. = 0. Hence from (17) 


2 Bie ? 1 o* 2) Fv na . liam a 
2. + : ih tanh (¢ = da(¢*) + | tan (: Te dQX(¢*) 


a T. \ 


bho 


3 pe ie 1/2 
+ = | tanh™’ ( + 2) d0(o*) = 0. (18) 
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Equations (2), (17) and (18) yield 





U ~ B. +m mW) (f + Pye r a)" + a) d0(¢*) 
? ‘J (: x =e)" + a) ane»)} + O(1/| w |? 


Since there can be no circulation about the aerofoil it follows from this result that 


a—a ab\ . (# a ay" 
( _ * 
([ + | \e + af q) 4006") 
* a— ey" d2(o*) = 0 19 
+f @ +0(2=%)" aaey =0. a9 
This equation can be regarded as a closure condition (see the closure conditions given 


in [2]). 
On y = 0 Ea. (16) yields 


o/ —< 2 > - , ; ‘ - o* — a) - > ah a) 
Bee (/ 7 [) so ‘6 + alle — + ane") 


on” {= sie + 0) dX¢*), (20) 


+E (a + $*)(6 — a) 


T Jane 


when —e < @ < —a,a < ¢ < Bb, and 


6(d) i. — : ( [ | + ') tan’ ie ° — oe + = d0(o*) 


* (¢* + a)(a — ¢) 
~azf tanh? ({G— 9G + o . 
R- tanh (e- + ¢%a— a} dQX¢*), (21) 


when —a < ¢ < a. (‘®’ denotes the ‘real part of’.) 

Two further auxiliary equations of less importance can be derived from the condition 
that on a well-designed aerofoil there should be no adverse infinite pressure gradients, 
except at the trailing edge. This requires that dg/ds > — ©, or since from (1) and (2) 


dq _ dq OQ0G _ q dQ 
ds dQ 0d OS - B 0’ 
that 
dQ 
G 29 
pm < oo, (22) 


Such infinite pressure gradients can occur at points separating the two types of boundary 
condition, i.e. at points C and D of Fig. 1. From (20), when | ¢ | > a, 


dQ I ib — a\” ( \(e 4+ ae Pas aa" 
0d - T(\p — QA) 2 + 2 { { . + / o* —$¢ o* + a d6($*) 


"(4 +2)(2= #1)" aan} 
” ee a+ ¢* amo") f, 


o*=-a 
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so that if (22) is to hold when ¢ = a, and when ¢ = —a, we must have 
( _+ / be + a) da(g*) — I (i=) ae") <0, (23) 
and 


(0+ NGA) een + (GS) aan =o, an 


respectively. This completes the mathematical theory. 
3. Method of designing an aerofoil. For a rounded-nose aerofoil the discontinuity 

in 6 at the leading edge, ¢ = —e, is 4x, while at the trailing edge, ¢ = b, it is 7, where 

7 is thus the trailing-edge angle. If for simplicity we assume that no other discontinuities 

in 6 exist, then in —e < @ < —a, anda < ¢ < b, we can write for the value of d@(¢*) 

occurring in Eqs. (18)-(21) 

dé 08 do* (25) 


. d6(o*) = d¢* — D ’ at 
ig i q 


~ ds Op 
where RF is the radius of curvature of the surface. 

Equation (21) can be regarded as the ‘design equation’ as it determines the aerofoil 
shape in the interval where the pressures are specified. However before it can be used for 
this purpose, it is clear from (25) that it is necessary to know g(¢) in —e < @ < —a, 
and a < ¢ < b, and this can only be found by first solving Eq. (20). It is apparent 
from (9) and (25) that (20) is an integral equation for qg. It can be solved by iteration as 
follows. We assume some likely distribution ¢(@); use it in 


i“ (26) 


to find s(¢), and hence from the known R(s), to find R(@). These values of R(@) and 
q(¢) are then used in (25), (18) and (20) to deduce Q(@). The first iteration is completed 
by calculating a new q(¢) from (9). When these iterations have convergedf, the final 
values of R(¢) and g(¢) are used in (21) and (25) to complete, without further iteration, 
the design of the aerofoil. Equation (19) can be conveniently satisfied by varying one 
of the parameters e, a or b during the iterative process. Finally the aerofoil co-ordinates 
follow from 
‘> cos 6 do ? sin 6 do m 
pf S—— gf = | (27) 
d q d q 


4. Special cases. (a) Flow about a given symmetrical aerofoil. 
From Fig. 1 it is apparent that the solution of this problem is obtained by letting ‘a’ 


tend to zero. In this case Eqs. (18) and (19) yield 


Sy 


| — da(g*) = 0, | o* do(¢*) = 0. 


¢%=—e Jo*=—e 


Subtracting (18) from (16), taking the limit as ‘a’ tends to zero, and making some use 
of (28) we find 
' 1 f° ; 
fw) = - | log (6* — w) dA(¢*). 
© detaxs 
tExperience with the same type of integral equation leads the author to believe that practical con- 


vergence is almost certain. 
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These three equations are the basis of the author’s method [4] of calculating the flow 
about given symmetrical aerofoils. 

(b) Usual equation of aerofoil design. This is arrived at by putting e, b = a in the 
above theory. If we write ¢* = —a cos y*. Then (18) and (19) become (see [1] and [2]) 


aT 


| Q(y*) dy* = 0, / Q(y*) cos y* dy* = 0, 
while on the aerofoil surface Eq. (16) reduces to the design equation [2] 


ical l ] 
6(¢) = es | O49") cot sy - 9 ~- ey” + »} dy*. 
T Jo 2 2 


(c) Riabouchinsky flow. Riabouchinsky flow [5] may be defined as flow past two solid 
bodies between which there is a bubble or cavity at constant pressure. All the steady 
symmetrical Riabouchinsky flows can be calculated from equations (16), (18) and (19) 
by putting dQ(¢*) = 0. (Unsteady Riabouchinsky flow is discussed in British A.R.C. 
Current Paper No. 149.) 

(d) A simple family of ‘roof-top’ aerofoils. ‘Roof-top’ aerofoils are those for which 
the specified velocity distribution is composed of two linear segments. Such aerofoils 
can be designed by putting dQ(¢*) = k,, (-a < ¢ < h), = ko, (h S @ S a), where k, 
and k, are constants, in Eqs. (18)-(21). 

For the simple example which has k, = k, = 0, and @ constant in —e < @ < —a, 
a < ¢ < b, except for discontinuities of 44 and 37 at the leading and trailing edges 
respectively, these equations yield 


m= —tant (F5) - Zente (SS7) 2 
_* (5 =)", (29) 
6) = x — tan PN — Lea PS aah 7 1 < bs OO 
and 
Q(¢) = 2, + R tanh {le + ie + Dy 2 
+ 76 tanh” (eve + by l¢|>1, (31) 


where ‘a’ has been put equal to unity without loss of generality. With the aid of (29) 
we find that the left hand sides of (23) and (24) reduce to 1/(e + 1) + 1/(b — 1), and 


1/(e — 1) + 1/(6 + 1), respectively. Thus, since e, b > 1, it follows from (23) and 
(24) that while an infinite adverse pressure gradient does not occur at ¢ = —1, it is 
unavoidable at @ = 1. 

Typical values of r and b for a modern aerofoil would be r = 12°, and b = 2.5, and 
from (29) and (28), e = 1.0117 and 2, = —0.1287. Suppose M. = 0.7, then from Eqs. 


(9) and (10) we find that the constant velocity of the roof-top is g = 1.22U. We shall 
omit further details of this example, but it should be noted that the flat nose occupies 
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less than 0.4% of the total potential difference over the chord, so that it could be rounded- 
off with but little effect on the velocity distribution. Eqs. (28)-(31) thus provide a com- 
paratively quick method of designing roof-top aerofoils. 

An alternative method of design would be to let e tend to ‘a’ and then specify the 
velocity distribution in the interval between the nose and the beginning of the flat 
section of the velocity distributions so as to obtain the desired nose-radius of curvature [2]. 
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TEARING AND INTERCONNECTING AS A FORM OF TRANSFORMATION* 


BY 
GABRIEL KRON 
General Electric Co., Schenectady, N. Y. 


1. Introduction. The mathematical problem under consideration owes its existence 
to the increasing complexity of engineering structures, and to a need for the piecewise 
formulation and solution of very complex physical systems. The author believes that he 
has solved his mathematical problem by experiment (in the manner of Heaviside) at 
least for linear systems and has also solved part of his program for certain types of 
non-linear dynamical systems. He is, however, an electrical engineer and not a math- 
ematician and now wishes to submit his thesis to a mathematical audience in the hope 
of arousing interest for a search into the mathematical legitimacy of his reasoning. The 
method proposed for analyzing and solving large-scale engineering structures in easy 
stages actually works in practice, as many textbooks and articles by scores of independent 
workers—both in this country and abroad—testify. There is no doubt that arguments 
and discussions, leading to a more rigorous mathematical proof, cannot but help the 
extension of the method into the piecewise solution of large-scale non-linear problems 
as well. 

2. The practical problem. The practical engineering problem the author attacks 
is as follows. Given a complex physical system containing mechanical, elastic, electrical, 
thermal, etc., subsystems, how can the “equations of state’’—as well as the “equations 
of solution’’—or ‘‘formulas of solutions’’—be established in easy stages, without manipu- 
lating or even writing down the simultaneous equations for the entire system? The 
procedure may be summarized as follows: 

(i) Tear apart the given physical system into a convenient number of subsystems, 
possessing no material contacts or other linkages with each other. (This last considera- 
tion is only a convenience, not an absolute necessity.) The two terminals of any point 
of tear need not be stationary but may have translations or rotations with respect to 
each other. 

(ii) Establish and solve the equations of each subsystem separately, as though the 
other subsystems did not exist. (One may use any established method of solution for the 
subsystems; or one may further subdivide each subsystem and use the present method 
to find its solution.) The solutions may be in a numerical or analytical form. 

(iii) Then interconnect the equations of solution (or the equations of state) of each 
subsystem by a routine procedure, to arrive at the solution (or at the equation of state) 
of the original given system. 

(iv) The remaining work consists of solving for (or eliminating) the comparatively 
small number of constraint forces that appear at the cuts. 

3. The theoretical problem. Instead of starting with one given complex system, it 
is possible to state the general problem in a slightly different manner. 

Let it be assumed that the equations of solution of a large number of isolated systems 
are available. The problem is how to utilize the already available partial solutions, in 
order to build up the equat’ ns of solution of all possible supersystems that may be 
constructed by mere interes ction of the given subsystems. 


L 
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The following third formulation of the problem will be used presently for a physical 
proof. Assume that the equations of solution (or equations of state) of one of the super- 
systems have already been established. The problem is how to establish, in a routine 
manner, the equations of solution of any one of the other possible supersystems (built 
out of the same set of subsystems) by utilizing the equations already derived. 

In the proof given below, each subsystem will consist of one electrical coil and each 
supersystem will consist of one possible interconnection of five coils. In practical problems 
there is, of course, no limitation on the size, physical nature, or complexity of each 
subsystem. 

4. An invariant procedure. The proposal to tear a system apart into several sub- 
divisions and to try to interconnect their equations has occured to many workers in 
these fields. The novelty of the present method lies in the proposed systematic and 
general procedure for accomplishing the interconnection. In particular, the proposed 
procedure enables one to utilize the entire apparatus of tensor analysis to organize the 
large variety of physical concepts and that of matrix algebra to organize the large number 
of constants into manageable forms. 

At the very beginning of his researches into the above problem, the author realized 
that to accomplish his objective he must follow a hitherto untrodden path. In particular, 
he conceived the idea of treating the processes of tearing apart a physical system and 
interconnecting again the component parts, as elements of a “group of transformations” 
Cz, (to.be discussed later). Thus he proposed to formulate or solve a given complex 
physical system by the following tensor (invariant) reasoning: 

(i) Tear apart the given complex system into a convenient number of independent 
subsystems. 

(ii) Set up the “equations of state’ (or the ‘equations of solution’’) of each sub- 
system in the form of a tensor equation. In the case of a dynamical system (for example 
a group of rotating electrical machines) the tensor equation of state of each subsystem 
(each rotating machine) assumes the form 

di® i dx” (1) 


> = -B — — 
Ca R ast + Las dt + ae dt ’ 


where I',s,, represents a general asymmetrical affine connection, not related to the 
metric tensor L,, . For each subsystem the geometric objects and tensors of R.s , Lag 
and I,,,,. are independently calculated. 

(iii) Represent the interconnection of the various subsystems by a “matrix of trans- 
formation” C%. . This matrix is established by simply inspecting the relations of the 
variables 7* and i* that exist at the points of separation—before and after the tearing— 
asi" = C2, i*. The components of C’%. are often plus or minus unity in simple systems, 
but in rotating dynamical systems they are functions of space and time. 

The existence of such a non-singular C is the key to the method of tearing. The rest 
of the paper elaborates a physical line of reasoning to justify and construct this single 
concept for as large a variety of practical cases as possible. 

(iv) Transform each geometric object by the aid of C%. in a routine manner, following 
their laws of transformations. Thus R,-s , Las and Ts, , « of the interconnected 
system are established in a routine manner. 

(v) The equations of state (or equations of solution) of the resultant system assume 
in terms of geometric objects the same form as those of the component subdivisions. 
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For instance, the equations of state of the interconnected group of rotating electrical 
machines becomes analogously to Eq. (1) 





ie ed ty ee ‘. 2) 
ee eae re a (2 


5. Invariance of equations. It should be emphasized that one of the basic require- 
ments of the method of tearing is that the form of the equations should remain invariant 
in passing from each member of the torn-apart system to the resultant system. (A second 
basic requirement, namely the invariance of power, will be treated later on.) It must 
also be emphasized that the form of the equations also remains invariant if the component 
subsystems are interconnected in any other possible manner. Moreover the invariance is 
maintained, even if each of the subsystems is further torn apart into several pieces and 
each of these pieces is interconnected into any of a large variety of supersystems. 

The method of tearing assumes that between each supersystem and each set of sub- 
systems there exists a non-singular matrix of transformation C and that their totality 


forms a “group’’. That is: 

(i) Each C has an inverse. 

(ii) The product of any two C’s is also an element of the group. 

(iii) The unit element leaves a system unchanged. 

6. “Orthogonal’’ networks. Perhaps the most important concept in interconnecting 
piecewise solutions is the existence of the inverse of every matrix of transformation C. 
Because of this property it is possible to pass freely from the equations of state to the 
equations of solution and back again at any stage of the analysis. To show the existence 
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a) "Mesh" network b) “Junction” network jc) "Orthogonal" network 











Fia. 1. Three types of representation of the same network. 


of an inverse C for all engineering structures, the simple electrical network of Fig. la 


will be used as an illustration. 
The network contains five coils, each coil having an impedance Z (or admittance 
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Y). It is well known that the equation of state (or equation of solution) of this network 
may be written in two different manners: 
° rae ° “8 1s 
(i) By writing two mesh equations e, = Za? (Fig. 1a). 
(ii) By writing three junction-pair equations J* = Y“"E, (Fig. 1b). 
[2 [ E| 
e Z (3) 
I Y 
For purposes of tearing it will be assumed, however, that it is always possible to 
write for the network as many equations of state (or equations of solution) as the sum of 
the number of meshes and junction-pairs (that is, as many as there are coils). The 
additional large number of ‘orthogonal’ equations may assume either of the following 


two forms: 


E a? * 3 ol & 
(Orthogonality of the meshes and junction-pairs is proved in Ref. [1], page 974.) 
That is, every network and every dynamical system may be described in a more 
general manner by assuming not only the active forces (e or J) but also the external 
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Fig, 2. Six different interconnections of five coils (stationary, linear graph), 
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constraints (£ or 7) as variables and writing for the system as many equations of state 
(or solution) as there are total number of variables. In conventional analysis such a 
course is not followed, but in the method of tearing it is necessary to have at least some 
of the superfluous equations available. 

In order to have the extra equations and extra variables at hand, it is necessary in 
the method of tearing always to start the study with the original physical system or a 
model of it. The latter always imply all the above variables and all the above equations. 

7. Example of stationary networks. Perhaps the simplest possible illustration of a 
group of C is given in Fig. 2, in which five stationary, one-dimensional electrical coils are 
interconnected in six different manners. (The number of possible interconnections of the 
same five coils is much larger.) The coils may have asymmetrical mutual impedances 
between them, and the impedances may be non-linear functions of the currents. 

Let it be assumed that the orthogonal equations of state, or the orthogonal equations 
of solution, of one of the networks are given by Eq. (4) or by a non-linear generalization 
of it. The problem is to establish the analogous equations for any of the other networks 
by utilizing the already known equations (4). The author’s thesis is that it is sufficient 
to establish a non-singular matrix of transformation C between the known network and 
the unknown network. With C established, the rest of the work is merely a routine 
application of tensor methods. 

Selecting any two networks, for instance the two-mesh, one-piece network of Fig. 
2A and the three-mesh, two-piece network of Fig. 2F, the matrix of transformation C? 
existing between them may be established in the following way. 

Reproducing the two networks in Fig. 3, assume as many current-variables in each 
network as there are coils, namely five. That is, assume in Fig. 3a two mesh-currents 
(i* , 2’ ) and three junction-pair currents. (Of course, each set of variables may be selected 
in a variety of ways.) In Fig. 3b assume three mesh-currents (7°, 7” , 7° ) and two 
junction-pair currents. Plot the path of each current through the network as shown 
(Again the paths may be selected in a number of different ways.) 

Consider next each coil in succession and equate the currents flowing through it in 
both networks. For instance, considering coil Z, , in Fig. 3a the current is 7* and in Fig. 


3bitisi* — i —7i*’ — 7’. Equating the currents in each of the five coils in succession, 
i?" =i —7" —-iw' — i" 
. + 7° “4 i = ie 
Q — i” =i" (5) 
ae - i? 4 + i , 
e+e 4+ = ae al 


(The five coils are identified in Figs. 2A and 2F.) 

It should be noted that the left-hand relation may be looked upon as C4 7“, trans- 
forming the network of Fig. 2E (containing isolated coils only, the so called ‘‘primitive” 
network) into that of Fig. 2A. Similarly the right-hand relation may be looked upon as 
Cy i’, transforming the primitive network of Fig. 2E into that of Fig. 2F. That is, the 


above set of equations may be written as 


= CPi = Chi’. (6) 
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Fia. 3. Transformation between two networks 


From this the relation between i* and 7” follows as 


o“ = CzCht', (7) 
where C? = (C%)~*. (This inverse always exists.) Hence the relation between the two 
network currents of Fig. 3 is 

i“ = Cp’, (8) 
where 
Cr = (C4)"'C? , (9) 
/ 1 —1 —-1 -1 
—] —|] 
C4 = | =] l - eR (10) 
—1] -—] —|] 
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This is the non-singular matrix of transformation from the two-mesh, one-piece 
network of Fig. 2A to the three-mesh, two-piece network of Fig. 2F. With its aid (or 
with the aid of its inverse), any physical entity or equation (such as Eq. 4 or its non- 
linear generalization) related to Fig. 2A may be transformed in a routine manner to 
that of Fig. 2F, or vice-versa. 

8. Existence of the “‘primitive’’ network. When only one network is given (say Fig. 
2F) and it is necessary to establish its equation of state, or equation of solution, the 
author advocates that the given network first be torn apart into subdivisions for which 
equations of state can easily be established. These subdivisions happen to be the indi- 
vidual “coils” for which the equations of state may usually be obtained by mere inspec- 
tion. The collection of isolated coils is called the “primitive network” shown in Fig. 2E. 

The step from the equations of the primitive network, say e¢, = Zagi’, to those of the 
given network, e,- = Zag? , requires only the establishment of a C. Since in most 
practical problems the given network needs to be looked upon either as a pure mesh 
network or as a pure junction network, it is sufficient in such cases to establish a rec- 
tangular, thus singular C. Nevertheless the non-singular C always exists. 

Any argument or proof based upon the existence of only a singular C may corroborate 
the formulae of the author, but are unacceptable as proofs of his method. The similarity 
of the singular C used in transforming from the primitive network to a given “mesh” 
network, to the singular C arising in “going from branches to meshes’’, is purely coinci- 
dental. The branches do not possess an equation of state of the form e, = 2Z,si°, but 
something more complicated. On the other hand, the method of the author is restricted 
to the use of only one single invariant equation, under any type of physical or hypothetical 
transformation. 

9. More complicated examples. A more complicated group of C exists in the presence 
of rotating electrical machinery, Ref. [5]. The ultimate element into which a system can 
be torn apart now is a two-dimensional cylindrical layer of winding (Fig. 4), in which 
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F1q. 4. Two-dimensional] cylindrical windings of rotating electrical machinery. 


the current-density wave is sinusoidal in space along the cross-section. The currents are 
removed through stationary brushes or rotating slip-rings. (One element of 7* represents 
now an angular velocity w.) Fig. 5 shows three different interconnections of three such 
layers of windings (one stationary and two rotating layers). The simplest arrangement 








154 GABRIEL KRON (Vol. XIII, No. 2 


is shown in Fig. 5B (the “primitive” rotating machine) in which all reference axes are 
stationary, are at right angles in space, and all windings are electrically isolated. The 
author has used this machine as the starting point for the analysis and solution of all 
other rotating electrical machines used in industry. The equations of performance of 
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A) B) "Primitive" machine Cc) 























Fie. 5. Three different interconnections of three windings (rotating, two-dimensional graph). 


any machine may be found from those of the primitive machine by merely establishing 
a C between the two machines. 

In a stationary elastic beam structure (Fig. 6, see also Ref. [3]), the ultimate sub- 
division (one beam) has twelve degrees of freedom (six for each end). More complicated 
examples of engineering structures that could be torn apart into a collection of smaller 
structures, could be cited indefinitely. 

10. ‘Tearing’ as a sub-group of affine transformations. Of course, simultaneously 
with the above C—representing actual tearing—many other transformation matrices 
also have to be used, in order to introduce hypothetical and physical reference frames 
that facilitate the tearing and the analysis. A reference frame, in which the tearing 
assumes a simple form, is usually not the best reference frame for the analysis itself. 
Practical exigencies of a problem (like knowing only real power in an electrical trans- 
mission system, instead of current or voltage) force the introduction of still other sets 
of reference frames. 

Most of the mathematical endeavors of the author have been concentrated upon 
fitting many of these physical and hypothetical transformations into the framework of 
tensor calculus. The need for considering ‘tearing’ and “interconnecting”’ also as a 
subgroup of the group of affine transformations is thus dictated by the appearance of 
so many other types of more conventional affine transformations in the analysis of complex 
physical systems. 

11. Practical considerations. In order that the solution of supersystems be in a 
practicable form, it is also necessary that all partial solutions of the subsystems involve 
far fewer non-zero elements than do the conventional solutions. Furthermore, the 
method for interconnection of the partial solutions should be comparatively simple and 
fast. It is believed that the suggested procedure satisfies these auxiliary—nevertheless 
important—practical considerations. 

To avoid any misunderstanding, the word “solution” refers not to a particular 
numerical solution of the variables, but to a general solution in which the variables still 
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occur in a general unspecified form. For instance, E, = Zasl° (a matric equation) is 
understood to be a general solution of J* = Y**E, if the elements of Z., alone are 
numeric. For any particular value of J“ a corresponding set of E, may be calculated. 
That is, the method of tearing is essentially an inversion procedure and the method 


should be evaluated accordingly. 
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A) Three different interconnections of seven beams 
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B) Stationary, six-phase, linear graph of a structure 








Fic. 6. Elastic structures. 


It will be shown in another publication that the analytical solution of problems, with 
all parameters factored out, also becomes feasible by the use of the method of tearing. 

In addition it must be emphasized that it is not the intention of the new method to 
compete with established methods of inversion (much less with methods of solution). 
The suggested procedure starts rather with “first-stage” solutions, already arrived at, 
and attempts to build up in successive stages large supersystems and their solutions. 
The new method intends to start where other methods end. Any overlapping in the first 
stage of solution is purely an educational or developmental accident. In presenting the 
subject this overlapping is unavoidable. 

12. Physical versus geometrical model. When only a set of equations is available 
(say a set of partial difference equations, or the dynamical equations of Lagrange) then— 
in order to use the method of tearing for solving the equations—the author’s first step 
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is to replace the set of equations by a physical model and then tear the latter apart, 
instead of the equations themselves. The author prefers to replace the given equations 
by a linear or higher dimensional graph (topological model) and to use electrical-circuit 
terminology in tearing and analyzing the model. Of course, the geometrical model (if 
correctly constructed) also contains all the superfluous variables and equations that the 
corresponding physical system does. The electrical-circuit exists only on paper and in 
general cannot be (and need not be) physically realizable. 

The electrical-circuit appearance and terminology of the author’s models should not 
mislead the mathematician into believing that he must become an electrical engineer in 
order to follow the author’s reasoning. On the contrary, it is the author who tries to 
become a topologist and solve topological problems of ‘‘cells’’ and ‘‘chains’’ as best he 
can. The author believes that the essentials of his method could be expressed in the 
language of linear and higher dimensional graphs (in the manner of Weyl, Ref. [10]) 
without any electrical connotations.* 

13. Tearing versus partitioning. It should be especially noted that it is not the set 
of equations that is being torn apart, but the physical system itself, or a model of it. 
Consequently the method of tearing has nothing whatever to do with the “partitioning” 
of matrices, which tears the equations themselves apart and not the physical system. 
The mutually exclusive nature of the two methods is clearly illuminated by the fact 
that the partitioning of matrices does not contribute any new information about the 
physical system, but the tearing apart of a model does. The method of tearing introduces 
additional variables,—namely the constraints (forces and velocities) appearing at the 
points of tearing—and their corresponding equations. (A more detailed comparison 
appears in Ref. [3]). 

It is the appearance of apparently superfluous variables and superfluous equations, 
that differentiates the method of tearing from other conventional methods of solving 
physical problems. Surprisingly enough, these extra variables and equations do not 
complicate, but rather simplify, the solution of physical systems and to a far greater 
extent than one would expect. 

14. Non-linear systems. For many years the author tried to achieve the inter- 
connection of the “equations of state’”’ of subsystems [5-9]. It is only during the last 
few years that he has also experimented systematically with interconnecting the “‘equa- 
tions of solution”’ of subsystems. In interconnecting equations of state no basic difficulty 
existed in dealing with either linear or with certain special non-linear dynamical systems. 
Similarly, in interconnecting equations of solution, no basic difficulties arise with linear 
dynamical systems. 

The author has published only one simple example for the interconnection of piece- 
wise solution of non-linear systems [4]. This example involved the inversion of Taylor’s 
series with several variables and its interconnection. Solutions for elastic beam structures 
(such as Fig. 6) with non-linear (plastic) characteristics in each beam have also been 
obtained by successive approximations. However, many difficulties lie ahead in inter- 


*Since the writing of this article the Author’s attention was called to a paper by J. P. Roth entitled 
“An application of algebraic topology to numerical analysis II. The validity of Kron’s method of tearing’’ 
(to appear in the Proceedings of the National Academy of Sciences). The paper identifies the following 
concepts of combinatorial topology that connect with the method of tearing: “Homomorphisms of 
homology and cohomology sequences induced by simplicial mappings of 1-dimensional complexes.” 
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connecting the solutions of more general types of non-linear systems. It is hoped that if 
competent mathematicians shed more light upon the pathway the author has been 
treading, the interconnection of many types of non-linear solutions will be facilitated. 

15. The invariance of power. The requirement that the form of the equations of 
each subsystem and of each interconnected system be invariant, is not sufficient to 
establish the laws of transformation of the various geometric objects. An additional 
requirement is needed. The author has made the discovery that a linear form, namely 
the total power input 7*e, (the product of generalized forces e, and generalized velocities 
i“) also remains invariant when a given number of subsystems are interconnected into 
any possible super-system. Because of the existence of the relation i*e, = i* e,- between 
any two systems, the usual laws of transformation of all tensors and geometric objects 
follow automatically. 

The above invariance is another key to the theory of tearing. The relation must 
exist, since the method of tearing does work in practice. But to prove the invariance of 
power in torn-apart systems in an unequivocal, scientific manner, acceptable to math- 
ematicians, is not a simple matter. 

The author has a plausible proof that appeared in Ref. [4], pp. 412-414 and 428-435, 
valid for stationary electrical networks only. If the proof is correct, a generalization of 
this proof to general physical systems should follow without difficulty. Naturally the 
proof is only a physical proof and not a mathematical one. Since the author is not a 
mathematician, he has never felt qualified to give a mathematical proof. However, all 
the assumptions he has made to prove his method to his own satisfaction have brought 
additional, useful contributions and so he feels that his proof is at least along the right 
track. In the following the author’s physical proof will be outlined. 

16. More general equations. In order to prove the invariance of power in any 
dynamical system, it is assumed that the most general form of the equations (of state 
and solution) must contain sets both of impressed and constraint forces (e, #) and 
velocities (7, 7) along each degree of freedom. For instance, for Fig. 7, the equations given 
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Fic. 7. Three points of view of a general network. 
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in Eq. (4) become generalized to: 


i+, +1.) | lr +E, e + £,] 
é+EH#, 2 Zo i+], Y, Y; (11) 
eo + E, 23 &4 12 + I; Y; Y, 


If the sets of mesh and junction equations are combined into one set, the above 
equations of state and solution may be written in the form 


Ca E. = Zas(t Pr. 
+ a(t + (12) 
;* 4 I Fo y7* E; + es) 


containing all possible forces (e, and J*) and constraints (7* and F,) along each of the 
topologically possible reference axes. In a similar manner, the dynamical equation (1) 
may assume the most general form 
8 8 
en tE, = Rat? + 1°) + Las ar A+ rT, (f% + PG" +1), (13) 
where a and 6 denote not only the customary generalized axes of active forces, but also 
their dual axes of constrained forces. 

Of course, it is never necessary to write down all the above equations with all the 
variables. Constructing a correct physical model is equivalent to having available (having 
written down) all the above equations in a tensor form. During the process of tearing 
and interconnecting the model, the engineer is enabled to write down only the absolute 
minimum number of equations and variables that he happens to need and to write them 
down only when needed. 

17. Invariance of total power. In the presence of the four sets of variables, it is 
always possible to view every orthogonal network (Fig. 7a) in two other ways: 

(i) As an “all-mesh” network (Fig. 7b) in which every impedance Z is short-circuited 
upon itself through a voltage, and forms a mesh. 

(ii) As an “all-junction” network (Fig. 7c) in which every admittance Y is open- 
circuited and forms a junction-pair. 

No matter what configurations are formed from the same n coils, the total current 
i* + J* in each coil always maintains itself constant (since each coil always remains 
shorted through the same voltage). Simultaneously, the total voltage drop e, + EF, 
across each coil also always maintains itself constant (since each coil always remains 
open). Consequently the total power in each coil—and in the entire network—remains 
constant, no matter what networks are built out of the same coils. That is, for all possible 
configurations of the same n coils the following relation is satisfied 


(ce, + E,)(i* + I*) = (.. + E, i” + 1°). (14) 


It should be recalled that the non-singular mairix of transformation C between two 
networks was established in Eq. (5) by assuming that in each coil of both networks the 
total currents are equal. 

This invariance of a linear form and the invariance requirements of the form of 
Eq. (12), enable one to establish the law of transformation of all tensors and geometric 
objects that arise in the method of tearing. Thus the author’s thesis, that the processes 
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of tearing and interconnecting physical systems can be viewed as a subgroup of the 
group of affine transformations, appears plausible. 

From a practical point of view the above thesis not only opens up the possibility of 
the piecewise analysis and solution of very large and otherwise unmanageable engineering 
structures, but it also allows one to concentrate all the resources of the calculus of tensors 
upon such studies. 

18. Acknowledgment. The author wishes to acknowledge his indebtedness to 
Professor Banesh Hoffmann for his constant encouragement and for many conversations, 
arguments and suggestions throughout the years, concerning the author’s endeavor to 
apply highly abstract tensor concepts to the analysis and solution of practical engineering 


problems. 
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Aerodynamics: selected topics in the light of their historical development. By Theodore 
von Karman. Cornell University Press, Ithaca, N.Y., 1954. ix + 203 pp. $4.75. 
When the leading authority in any field writes a book in that field, one looks forward with anticipa- 

tion to the reading of that book. The book under review here is just such a volume. In fact, its history in 

my Office will probably be many times repeated. One day it arrived and lay unwrapped on the secretary’s 

desk. Each person passing her desk stopped, looked and then picked up the book and sat down for a 

few minutes leafing through to get an idea of what was covered and how it was covered. In the present 

case, the book then disappeared, to turn up a day or two later in an adjacent office where it had been 
taken for more careful examination. However, I as the intended reviewer, got hold of it and read in detail 
the things that others were attempting to gather by hasty examination. 

The book is not a text. It is intended as a presentation of the physical aspects of the ideas of aero- 
dynamics and the historical order of development of these ideas with a few pertinent remarks about the 
men who made those developments. The very early aerodynamic ideas before the period of flight are 
shown in a perspective only possible now that the correct ideas can be selected from among the early 
incorrect ones. It is pointed out that some of the basic ideas are very early indeed. The development of 
aerodynamics ideas since actual flight of man began is traced through the lifting line notions, as developed 
by Lanchester and Prandtl. The importance of the ideas of Kutta and Joukowski are noted, and the 
general consequences and verification of the circulation picture as it developed is presented. Finally, 
some of the very recent ideas on theory of low aspect ratio wings as developed by Jones are described. 

Next, the long, slow development of ideas about drag and skin friction which involves the develop- 
ment of the ideas of boundary layer, the vortex streets and the relation between laminar and turbulent 
flow are all given in relation to the contemporary flight problems. Next, a section on supersonic aerody- 
namics starts with the usual elementary picture of moving sound sources and goes on to develop linearized 
wing theory to introduce the shock wave. All these phenomena together with shock wave boundary layer 
interaction and a development of swept-back wings are given in relation to modern high speed flight. 

In Chapter 5, a section on stability and aero-elasticity presents the dynamics of flight and the related 
problems of the stability of the structure of the aircraft. It is pointed out how the various instabilities of 
the flight path and instabilities of the elastic structure were gradually learned and understood and cor- 
rected. Finally, the aerodynamic problems of propulsion are considered in a chapter with a title suitable 
to the book at hand, ‘‘From Propeller to Space Rocket.”’ As one reads this last section, one is impressed 
with the similarity between our understanding of space travel and the understanding of air travel as of, 
say, 100 years ago. No doubt, nearly all of the necessary basic ideas for space travel are already included 
in the many ideas that have been treated or speculated upon, but it remains for the future to show us 
which are the good ideas and which are not so useful. 

In the propulsion section, general propeller theory and the jet and rocket engine types of propulsion 
are presented so far as basic ideas are concerned, and are given with proper reference to the originators 
and those who first worked on the various types. 

Throughout, the book is written to appeal not only to professional aerodynamicists, but also to any- 
one who wants a good physical picture of where things stand in aerodynamics, and how we got there. 
Throughout, the author has introduced not only historical information and serious discussion, but also 
the occasional witty remark or wisecrack which has become popular in the aerodynamics field or is worthy 
of that fate. 

The book will be of value both to the aerodynamics expert and to the educated layman and both will 


find it well worth the time of reading. 
Howarp W. Emmons 


Introduction to elliptic functions with applications. By F. Bowman. John Wiley & Sons, 

Inc., New York, 1953. 115 pp. $2.50. 

The purpose of this book is to present a short, practical, elementary account of the properties of the 
Jacobian elliptic functions and their applications for physicists, engineers and applied mathematicians. 
The book, unfortunately, is written in the vein of an elementary undergraduate text and includes so 

(Continued on p. 168) 
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TEMPERATURE DISTRIBUTION IN SOLIDS OF VARIABLE THERMAL 
PROPERTIES HEATED BY MOVING HEAT SOURCES* 


BY 
R. J. GROSH, E. A. TRABANT, anv G. A. HAWKINS 


Purdue University 


Abstract. The temperature distribution in rods, plates, and infinite solids heated by 
moving plane, line, and point heat sources is obtained for the quasi-stationary state. 
The thermal properties of the solid can vary as a linear, parabolic, or some other con- 
venient function of temperature as long as the thermal diffusivity of the material remains 
nearly constant. 

Introduction. The theory of moving heat sources has been applied in the main to 
welding problems; however, in recent years the theory has been extended in applications 
to extrusion, flame hardening, gun barrels, and determinations of thermal diffusivities. 
The theory as usually presented assumes that thermal properties are independent of 
temperature, that heat generation due to phase change or Joule heating is negligible, 
and that the source moves in such a manner that the temperature field about the moving 
heat source is independent of time, the so-called quasi-stationary state [1].' It is the 
purpose of this paper to present a theory which accounts, at least in part, for the varia- 
tion of thermal properties with temperature. The results obtained are applicable to 
certain one, two, and three dimensional problems where the heat source can conveniently 
be represented by a moving plane, line, or point. 

Fundamental equations. The partial differential equation for heat conduction in a 
moving isotropic solid the thermal properties of which are functions of temperature 
and/or position and where any heat generation within the solid is negligible can be 
written as 


oe DT 
V-(kVT) = pe, (1) 
Dr 
where 7’ is the temperature and a function of the coordinates 2, y, z in a fixed Cartesian 
coordinate system and time 7; k is thermal conductivity and a function of temperature 
and/or coordinates; p is density and a function of temperature and/or coordinates; 
c is specific heat and a function of temperature and/or coordinates; and DT’/Dr is the 
substantial, fluid, or total derivative of temperature. Assuming linear motion for the 
solid as shown in Fig. 1 and a quasi-stationary temperature field, this equation can be 
reduced to 
oT 
V (kVT) = —pev—, (2) 
Ox 
where v is the speed of the solid with respect to the fixed Cartesian coordinate system. 
The source will be considered fixed as a plane, line, or point about the origin of the fixed 
coordinate system. It is of no consequence whether the source is considered fixed in 
space and the solid moves, or, the solid is considered fixed in space and the source moves. 
*Received July 19, 1954. 
1Numbere in brackets refer to the Bibliography at the end of the paper. 
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a 


In the latter case, the fundamental equation may be expressed as Eq. (1) through the 
use of a suitable transformation of variable [3]. 

Heat sources are assumed to be planes, lines, or points. A plane source is a plane 
surface of area A which emits q, units of energy per unit time in its normal direction n as 


q, = —kA ee (3) 


2.38mm Diameter Welding Rod 
Direct Current Arc: 85a. ,19.5v 


Stainless Steel Plate 





™ Pee ee 
Fixed Cartesian Re(x'e yo+z") 


Coordinate System | 








Solid Assumed 
Moving in This Direction 





Fic. 1. Coordinate system. 


A line source is a segment parallel to the z axis of length L which emits q; units of energy 
per unit time as 


: oT 
qx = Lim — 2aLrk —, (4) 
r-0 or 
where r = (x° + y’)'””. A point source is defined to emit g, units of energy per unit 
time as 
. rk - 
gq, = Lim — 4rh*k —, (5) 
R-0 Ok 


9 > V1, 


where R = (& + y + 2) 
Other boundary or surface conditions will be considered to be of the adiabatic or 
“no flux” type. That is 


> 


oT : = 
aE = ( Eg=+0, +b, (6) 
which follows from 
9T a 
k—=0 £=+0, 4b, (7) 
0& 


since k > 0. The symbol é represents z, y, z, r, or R; b is a constant. 
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The thermal properties of a solid are generally functions of temperature. It is con- 
venient to assume that for those properties appearing in Eqs. (1) through (7) 


k = kf'(T), (8) pe = (pc)of"(T), (9) 


where k, and (pc), are positive constants and f’(7) is the temperature derivative of 
some function of temperature f(7). It is assumed that f(7) and f’(T) are differentiable 
and f’(7’) is greater than zero to avoid violation of the Second Law of Thermodynamics. 
The particular convenience of these forms appeals, in any case, to experimental data. 
The efficacy of this form of assumption is shown in Fig. 2 for a stainless steel. 
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Fic. 2. Thermal conductivity, product of density and specific heat, and thermal diffusivity. © Data of 
Powell [2], CGS units. 


Generally, 


ee 
. * f(T) 3E" (10) 
This equality is used to simplify Eqs. (2) through (7); they become respectively 

V'° S(T) = —2nrv of(T) (11) gq. = —kA of) (12) 

Ox on 
gq, = Lim — 2ark,L np (13) 

y o>, Of(T) 
, = Lim — 42R*k, , 14 
eS — ae ied 
af(T) s) _ 

._ > 0, £= +0, +b, (15) 


where \ = 3a and a is the thermal diffusivity k/pc. 
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The reader is referred to papers by Rosenthal [3] and Rosenthal and Cameron [4] 
for solutions of the following boundary value problems after the above transformations 


have been made. 
One dimensional heat flow problem. If it is assumed that heat flows in the x direc- 


tion only, then 


0 earn 23) —— oT f 
Ox {he (7) Ox) a —(pe)ot J (7) ox (16) 
Let the boundary conditions be 

aT 
ite =W6 oman. (17) 

Ox 

oT 
—kA— = 4q, = 0, 8) 
ax q x (18) 


and let the initial temperature of the solid be 7’ : 


T = T; r= +o, (19) 
From the foregoing transformations, this problem reduces to 
f(T) _ _o, 2f(P) 
ax" — oe. 
of(T) 
“a =0 «= +2, 
Or 
(20) 
of(T) 
—kA—— = ¢ r=0 
Ox 
{(T) = f(T,) r=+a 
The solution to this problem is 
ae — q: ' 
{((T) — f(T; = ‘exp (—2dvz) gz > ®, (21) 
, A(pc) ov 
es ray er / 7] ¢ a 
{(T) — f(T) = 4+ x <0. (22) 
. A(pc) 1 — 
If f’(T) = 1 + mT, where m is a positive constant; that is, the thermal conductivity 


and product of density and specific heat are linear functions of temperature; then Eqs. 
(21) and (22) reduce to 
7 ] f 2q m , m \2 i ) P 
r =—% exp (—2dox) + (1 + mT;) —1 r>0, (23) 
m | A(pe)ov = 


nr 1 ff 2am ,, ae 
r =i Oh ++ mr | — 1 z= 0. (24) 
m \L A(pc)ov f 
The solution corresponds to the temperature distribution in an infinite rod or solid 
due to a moving plane source. 
Two dimensional heat flow problems. (a) The infinite plate. If heat flows in the 
xz and y directions from a line source along the z axis, then 
a f(T a’ f(T af(T) 
IAT) , FI) y ———. 


“= —2Xr 25) 
Ox” oy Ox 38) 
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Further, let 


of(T) 
oe] a 
or 0 ‘ ‘ (26) 
gq. = Lim — 2arLky ad (27) 


and the initial temperature of the plate be 7; so that 


(T) = f(T,) r= o, (28) 


The solution to this problem is 
f(T) — f(T) =; = exp (—Avx)K,(dvr), (29) 
: QrkoL 


where K,(dvr) is a modified Bessel function of the second kind, zero order, with argument 
dvr. If the thermal properties are certain linear functions of temperature; that is, 


f'(T) = 1 + mT; then 


T = ‘| 2 exp (—dvx)K,(dor) + (1 + mr | "i i (30) 
m alk, 
1200F 
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Fic. 3. Application of theory Eq. 30 to experiment. 


This equation is of interest in regard to welding thin, stainless steel plates. An applica- 
tion of this equation to welding can be found elsewhere [5]. A synopsis of this is shown in 
Figs. 1 and 3. 
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(b) The finite plate. In contradistinction to the previous case, if the plate is of width 
b in the y direction then 


of(T) 
= es i) = 0, db, ‘ 
oy y ‘ (31) 
f(T) , 
oJ\!) = () r= +o. (32) 
Ox 
The solution to this case, obtained by Rosenthal using the method of images, is 
(TT) — f(T, = i. exp (—)Av2) 2. K,(rur,), (33) 
where 7, = (x + [y + 2nb]’)'”’. Again, considering the linear type variation of thermal 
properties, 
5 1 )} 2qzm = PY us 
T= — 4 — exp (—)vz) K,(dor) (1 + mT’;)° — 1?. (i 
= | 2a I 2 { (Avr) + + mT;) 1 34) 


Three dimensional heat flow problems. (a) The infinite solid. If heat flows in the 
x, y, and z directions from a point source moving along the z axis, then 


V°s(T) = —2rv of) (35) 
: Ox 
The boundary conditions will be 
“oe =0 R= (36) 
ge = Lim — 4eR’h, ae (37) 
and the initial temperature of the solid will be 7’; or 
f(T) = f(T) R= o, (38) 
The solution to this problem is 
f(T) — f(T) = dr exp (—Avzx) exp (—)vR). (39) 
4rk,R 


If the thermal properties are certain linear functions of temperature, then 


1 ial f na 1/2 
T = - 4 EF exp (—dAvax) exp (—AXR) + (1 + mr) | — vs (40) 


(b) The finite solid. If the solid is of infinite extent in the x and y directions, but 
bounded by the planes z = 0 and z 


g; that is, the boundary conditions are 


of(T) 
an ee Tr co, (41) 
or 
af(T 
os f= 9 s=0, 9 (42) 
Oz ‘ 
K(T) = f(T) r= © 


(43) 
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then it can be shown that the solution to this problem is 


f(T) — (7) = — exp (—rz) © exp (NE) (44) 


2rkp 


where R, = [x? + y* + (2 + 2ng)’]'”. 

(c) The solid bounded by a cylinder. If the solid is of infinite length in the z direction 
and is bounded in the y and z directions by a cylinder of radius a concentric with the 
x axis, another solution is available in rather complicated form. The cylindrical surface 
is assumed insulated and a point source moves along the z axis. 


Then 
V'°f(T) = —2dv on (45) 
of(T) = () g _ (y? +- g*)'4 =a, (46) 
0& 
af(T) 
of(T) =0 r= +0, (47) 
Ox 
aiid 2p, f(T) 
q, = Lim 4rR°k aR” (48) 
f(T) = f(T, r= +o, (49) 


The reader is referred to the original paper of Rosenthal and Cameron [4] for a discussion 
of the appropriate solution. 

It is to be noted that, in all the above cases of one, two, and three dimensional con- 
duction, only linear types of properties were considered. The method outlined is able to 
account for any variation of properties as long as Eqs. (8) and (9) are valid. However, 
the use of more accurate approximating functions requires solutions of high degree 
algebraic or trigonometric equations. 

Summary. A method is outlined to account for certain variations of thermal proper- 
ties of a solid when heated by a moving heat source. Variations considered are of the 
type where thermal conductivity and the product of density and specific heat are similar 
functions of temperature; or, where these properties vary but their ratio, the thermal 
diffusivity, is constant. The theory is applicable to one, two, and three dimensional 
cases of heat flow from a plane, line, and point moving sources. 

It is felt that the theory furnishes a useful approximation to certain real systems. 
The theory has been verified in the particular case of welding thin stainless steel plates. 
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much detail and repetition that anyone competent in mathematical manipulation could more easily and 
quickly obtain the desired material from a standard handbook while the less sophisticated person would 
find it difficult to draw the needed information from the detailed mathematical computations. Most of 
the formulas in the book are developed in exhaustive detail and the same applies to the exceedingly 
numerous illustrative examples. There are a total of 128 problems for the reader to work out, two-thirds 
of which are relatively trivial exercises in manipulation. No numerical tables are given but references to 
available tables are made throughout the book. (See, however, A. Fletcher, Guide to Tables of Elliptic 
Functions, Math. Tables and other Aids to Computation, 3, 229-281, 1948.) Except for the above men- 
tioned shortcomings, this book does present as much information about Jacobian elliptic functions and 
integrals as the average consumer would ordinarily need. The collection and detailed analysis of conformal 
transformations in which various elliptic functions and integrals appear certainly illustrates the usefulness 
of these functions, in particular in relation to the use of the Schwarz-Christoffel transformation and the 
chapter on the reduction of elliptic integrals to standard form is clearly presented. (For a better and more 
useful presentation see, however, P. F. Byrd and M. D. Friedman, Handbook of Elliptic Integrals for 
Engineers and Physicists, Springer-Verlag, Berlin, 1954.) 


PETER CHIARULLI 


Proceedings of the symposium on nonlinear circuit analysis. Sponsored by the Microwave 
Research Institute of Polytechnic Institute of Brooklyn, in cooperation with The 
Institute of Radio Engineers Professional Group on Circuit Theory, and co-sponsored 
by the Office of Naval Research, The Office of Scientific Research, The Signal Corps, 
New York, 1953. xi + 411 pp. $4.00. 


This book, the second volume to appear in the Network Symposia Series sponsored by the Micro- 
wave Research Institute of Polytechnic Institute of Brooklyn, concerns nonlinear network analysis, i.e., 
nonlinear differential equations. Following several introductory papers by Weber, Stoker, Friedrichs, 
and Lefschetz, are seventeen papers covering topics from purely numerical analysis of nonlinear differ- 
ential equations to physical systems exhibiting subharmonic behavior. 

While a complete table of contents is unwarranted here, a breakdown of the classification of papers is 
in order. Excepting the four introductory papers, six treat mathematical techniques for analyzing non- 
linear systems, two each on perturbation methods and response of systems to special inputs, one each 
on stability of differential-difference equations, impossible behavior of nonlinear networks, inertial 
parameters, synthesis of nonlinear systems, subharmonics, magnetic amplifiers, and intentional non- 
linearization of servomechanisms. 

The format of the Proceedings is good and, while many minor errors were found, none of them would 
be a hindrance to the serious reader. This volume is a desirable addition for those working in nonlinear 
mechanics and is even recommended to those who would like to familiarize themselves with the field. 


SHELDON LEvy 


Transactions of the symposium on fluid mechanics and computing held at New York Uni- 
versity, April 23-24, 1953. Edited by G. Birkhoff, K. O. Friedrichs, and T. E. Sterne. 
Interscience Publishers, Inc., New York, 1954. iv + 243 pp. $5.00. 


The volume contains the papers presented at a symposium that was sponsored by the American 
Mathematical Society and the Office of Ordnance Research, U.S. Army. The contributors to the volume 
are M. J. Lighthill, G. F. Carrier, G. Birkhoff, G. S. S. Ludford and M. H. Martin, J. H. Giese, L. Bers, 
A. Weinstein, P. Germain, R. v. Mises, M. Lotkin, P. D. Lax, L. H. Thomas, G. E. Hudson, and 


H. Schardin. 
W. PRAGER 
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ON TRANSVERSE VIBRATIONS OF THIN, SHALLOW ELASTIC SHELLS* 


BY 
ERIC REISSNER 
Massachusetts Institute of Technology 


1. Introduction. This note is concerned with the form of the differential equations 
of the theory of shallow shells, as given by Marguerre [2], for the case of vibratory 
motion of the shell. In their general form these equations may be reduced to three 
simultaneous differential equations for three displacement components. It is our principal 
object here to show that for types of vibrations which we call transverse the equations 
of the theory—with negligible error—admit of an important simplification. This simpli- 
fication consists in the omission of longitudinal inertia terms. It is possible with this 
omission to treat transverse vibrations by means of the same system of two simultaneous 
equations for a stress function and for one displacement component which is known 
for the problems of statics of a shallow shell. 

The advantages of the present treatment of the problem of transverse vibrations 
are illustrated by determining the frequencies of transverse vibrations of a shell with 
middle surface given by a second-degree equation if the shell is simply supported along 
a boundary with rectangular projection on the base plane. 

We have previously used a special case of the results of this paper in the analysis of 
shallow arches [4]. Among other applications which are possible we mention the problem 
of the axisymmetrical vibrations of shallow spherical shells which has previously been 
considered by Federhofer [1] and the present writer [3], without taking advantage of 
the simplification which will be shown to be possible. 

2. Basic differential equations. Marguerre’s equations of the linear theory of 
shallow shells are here taken in the following form 


eee, a 
oat + 2 = oh * (1b) 
+ ees | 2 Ne +5 v..,| +z 2 New + > v,| = hoe, (2) 
aa + a — = 0, (3b) 
M, = - (2% +» o*), (4a) 
M, = —- D(a +y *v), (4b) 
Oy On 
- *Received August 5, 1954. A report on work supported by Office of Naval Research under Contract 


N5-ori-07834 with Massachusetts Institute of Technology. 
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°w 
M,, = —(1 —v)D — (4c) 
Ox OY 
Ou Oz Ow N, —vN, * 
= — —— = —— (5a) 
Ox Ox Ox Eh 
Ov Oz Ow N, — wN, . 
: Sa oa : , (5b) 
oy Oy Oy Eh 
Ou Ov Oz Ow ; Oz OW N., —— 
Oy + Ox 6 Ox Oy * dy dx Gh° ad 


In these equations 2, y are Cartesian coordinates, the equation of the middle surface of 


the shell being z = 2z(z, y); u, v, w denote the displacement components along the coordi- 
nate axes, ¢ the time, p the density, h the shell thickness, N the tangential stress resultants, 
Q the transverse stress resultants, M the stress couples, and D = Eh*/12(1 — v’) the 
bending stiffness. 

Introduction of Eqs. (3) to (5) into (1) and (2) leads to three simultaneous equations 
for u, v, w. Restricting attention to cases for which F and h are constant these equations 
may be brought into the following form 
i~-¢# Vu + 1+» 9 (2 4 a) 4 i= 2 (2 Tw + ow v%) 


2 Ox Ox 





2 2 x \Ox a 
ite sie (6a) 
es 1+vr 9 (2 ow dz 4 _ p Ou 
2 dx \dx dx dy dy/  ~E at’ 
1—v_p» ] v do fou v 1 — v (dz > DW —> 
Se Y= LE < [= + ) aaa — (: Vw + — V% 
2 2 ody \dz oy 2 oy oy (6b) 
(6b 
+ ee (= dw , a oe) _ pdv 
2 dy \dxdx ~ dy dy E ot” 
— Eh Xz | du dz Ow av dz dw 
—DV*V"w + {¢ a a | 2 5 ae ae 
1—v (\dx° Lox Ox OX \0y Oy OY! 
dz | av dz OW Ou , dz Ow 
+25! + > oe, (4 % au) 
Oy Loy Oy OY ON Ox Ox 4 
(6c) 
2 9 Av dz Ow , Oz Ow 
+a —9 fe [ey wy eae, oul 
Ox Oy LOY Ox Ox OY Oy Ox 


E dz Ou Oz a4 
= phi <3 —- = 33 — = 3 |- 
Ot Ox dl Oy oat 
Equations equivalent to these have been used by Federhofer [1] and the writer [3] for 
the problem of axisymmetrical vibrations of shallow spherical shells. It is an indication 
of the difficulties associated with this system of equations that no successful attempt is 
known of an evaluation of the frequency determinant obtained from this theory for 
what is probably the simplest of these spherical-shell problems, the problem of the 
complete shell segment with clamped edge. 

3. Simplified differential equations. It is well known that in problems of statics 
considerable practical advantage accrues if use is made of two simultaneous equations, 
for the displacement component w and for an Airy stress function F instead of three 
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simultaneous equations for three displacements. In problems of dynamics this approach 
is prevented by the presence of the longitudinal inertia terms u,, and v,, in Eqs. (1). 
Suppose, however, that these terms are negligible for certain classes of dynamic problems. 
Then use of Airy’s stress function becomes possible in the solution of these problems. 
Our principal aim in this note is concerned with just this point. We shall show that for 
classes of problems in which vibrations take place principally transversely we may 
neglect longitudinal inertia, at the same time when longitudinal straining may exert an 
important influence on the dynamic phenomenon. 
Omitting u,, and v,, in Eqs. (1) we set 
oF , OF OF 


N => > nn. = d y = — 
, Ney Ox Oy 


: —z, 7a, b, ec) 
oy” Ox” (@ 


and utilize a compatibility equation of the form 


oO (du Oz =) o (% Oz ~) a 2 ev (2 re dz Ow 
dx \dy © dy dy Ox Oy \Oy Ox Ox dy ~~ dy Ox (8) 


ov Oz Ow Oz oe) 
oy” \Ox Ox OX) 


2 2 2. «2 2 2 
Oz Ow Oz Ow 02z0W 


= 2 _-—s 7 > 
~ Ox Oy Ox Oy Ox Oy dy Ox 





Introduction of (5) and (7) into (8) gives the first of the two simultaneous equations 
for F and w, 


_ V2} =6Ow a2 Ow az Ow 
V' VF = El f, —— a ee ee oe oe eu I 
vv : ‘| Ox Oy Ox Oy Ox Oy” Oy Ox” () 


Introduction of (3), (4) and (7) into (2) gives the second of the two simultaneous equa- 


tions, 
— aw az OF 02 OF . a2 OF 
DV?V>w + ph S53 = —2 SS 4 4 GSS. II 
ils ot Ox Oy Ox OY Ox Oy Oy Ox (1) 





4. Justification of simplified differential equations. In order to establish conditions 
under which neglect of longitudinal inertia is justified we undertake an order-of-magnitude 
analysis of the displacement differential equations (6). We introduce dimensionless 


variables 


(9a, b, ce) 


wre 
I 
3 
I 
ik 
4 
€ 


and 


w= Wf, u = Ug, v = Vk, (10a, b, ec) 


such that the functions f, g, k and their derivatives with respect to £, 7, 7 are at most 
of order of magnitude unity. We further set 


H H 
Zr = L c. , zy = L ce ’ (lla, b) 
H H H 
@s2 = 73 Ces » 2w = 73 or Z2u = 7? Pita (12a, b, ec) 


— : ” ; . 
where ¢. , fy , fez » Sy» fx, are at most of order of magnitude unity and where (H/L) 


is negligibly small compared with unity, in order that the theory of shallow shells be 


applicable. 
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We now consider terms of all representative types in Eq. (6). From (6a), 


U H.W H W. pw 
= Gee a Fe fee + st hes a Se ooo  &—— Uog,, . (13: 
Je TT $z pie + 73 hs 1 Set E g 2) 
Equation (6b), being analogous in structure, need not be considered separately. Equa- 
tion (6c) becomes 


ee es Ce { U H.W a 
1201 — ) 7p | Jetee }+ > td] 5 J: + L Gs ] fe + 4 


(13b) 
= pha'| WS. _ . t.U9,, er | 


4 


We are free to dispose of four quantities in (13), U, W, 1 and w. Appropriate dis- 
position leads to the results which are desired. 

Since we expect that longitudinal straining will exert an influence on transverse 
vibrations we must see to it that in (13a) terms involving both f and g remain significant. 
This is accomplished by setting either 


U = L W (14a) 
or, if 1/L > 1, 
: By t.. 
U= LL W. (14b) 
This transforms (13a) into 
l : pw I” al 
Jee + Sefer + L fale Too = E 9" (15a) 
or 
f : iy . 
Oe + 2 babe + bach te = ae (15b) 


Since terms in (15) can be at most of order unity we conclude that (14a) and (15a) are 
applicable as long as 


; = (0(1), (including < 1), (16a) 
while (14b) and (15b) are applicable as long as 

L ; : ' ' 

F ie O(1), (including < 1). (16b) 


Case (16a) means that representative vibration wave lengths / are at most of the order 
of magnitude of representative shell shape wave lengths LZ. Case (16b) includes the 
case of vibration-wave lengths which are large compared to shell shape wave lengths. 
Case (16a) is of importance in problems such as that of the shallow spherical shell. 
Case (16b) concerns problems such as the problem of vibrations of corrugated plates. 
In either case it may happen that the coefficient on the right of (15), pw'l’/E, is small 
compared to unity. If this does happen then longitudinal inertia will be negligible. 
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To obtain information on this point we consider next (13b) with either (14a) or (14b). 
From (13b) and (14a) 


l + 1 Hh 
9 le tin Sen + *°*] + oo oar Selle + fatal + °° 7 
12(] Vv) l 1 v L l wy 
(17a) 
_ pha? fH ve] 
+ = E E ty CMe + ° 
From (13b) and (14b), 
l h° ‘ 1 HH’) L 
, 2) = [—Sece + °°*) + - ral a + $. Ls, | a ++} 
121—yr)l l-—y L l " 
(17b) 
= iP Poe 
+ — E E LF. L G29 rr 
Since transverse inertia is to be important, it is necessary to match the coefficient of 


f.. with at least one of the coefficients on the left side of (17). We consider first (17a). 
We set 


1 | pw = 
wi-r)t E (18a, i) 





or 


1 _ pw (18a, ii) 


l-y Ll E* 
With this, taking account of the fact that (H/L)* < 1, 
H’ I 


l~Seur + *°*} + 12 RE Secilge + Sefe) tees} tees = See (19a, i) 
or 

my — i” 

12 FE [—feeee + -°°) + Secllge + Sefe] + -°°} Hee = See (19a, ii) 


Equation (19a, i) is applicable as long as 
12 —3 +35 = O(1) (or < 1), (20a, i) 


while (19a, ii) is applicable when 


1 L* h? 


a ee i ° - 
12 PP O(1) (or < 1), (20a, ii) 


the two regions of parameter values complementing each other. 
We now introduce pw’/E from (18a) into (15a). The coefficient on the right of (15a) 


becomes 


1 h? : 
2d —v) F 
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and 
ae | . 
l—-ry LL (11) 


respectively. We may assume that we have, besides H’/L’ < 1 also h’/I’ « 1. In fact, it 
is only if this assumption is made that the effect of transverse shear deformation can be 
neglected, as we have done at the outset. Since //L = O(1), we have then that both 
the above longitudinal inertia coefficients are small compared with unity. Accordingly, 
as long as representative wave length | and frequency w are related as in (18a, i) or (18a, ii) 
and as long as 1/L = O(1) we have that the effect of longitudinal inertia in the theory of 





shallow shells is negligible compared with the effect of transverse inertia. We remark that 
further simplifications of the basic differential equations are possible if 

l Hy’ Fr 1 Dh’ 

7K1 and(or) 127555 « or 5 p<), 


provided, in the third case, that the boundary conditions are such that bending need 


not occur. 
We now consider case (14b), (16b) and (17b) for which //Z may be large compared 


with unity. We set in (17b) 


l h* pw ‘ 
= a (18 
21 -—-y) f° E 18b, 1) 
or 
l H* pw eal ts 
1 <r rs bait E . (18b, 11) 


This gives 


H? I Gs \ 
a ae Ome 9 -——— ; ba - me. 
[ J eeee ]+ 12 h- Py me + $z l fe| + J 
ee — 
+---=f,- ie Aaa tee (19b, 1) 
or 
[vi f . Bs 
12 H? 7* L—J eee + eee] + ral 0 + $z 1 fe| + i 


a. ” 
occ = —- ¢ — eee ( 
-+}- | ite zr G29rr ? (19b, il) 


while (15b) becomes 


l h? . 
I 12a — 7) PI ania 
Yee + fete + 7 Cfee tcc = 1 
id HT —— 
1 re y | ty Grr ° (21 ’ 11) 


Equations (19b, i) and (21, i) can be further simplified. Since 


ee. ee 


ML = wT = 2: 
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in the range L/] = O(1) and l/h > 1, it follows that (H/L)’(l/L) « 1 and the second 
term on the right of (19b, i) is negligible. Furthermore, the longitudinal inertia term 
in (21, i) is also negligible. Accordingly, the simplified equations (I) and (II) are appli- 
cable even when //L > 1, provided only that (18b, i) holds, that is provided that plate 
bending action contributes significantly. 

It is not possible to simplify Eqs. (19b, ii) and (21, ii) further since we are not allowed 
to conclude that the inertia terms (H/L)*(1/L)g,, and (H/L)’(l/L)*g,, are always neg- 
ligibly small. Accordingly, when 

+> —_— “3 a: x, = 0(1), 
the simplified equations (I) and (II) may not be applicable. 

5. Vibrations of paraboloidal shell with rectangular boundary. We consider as a 
simple illustration of the use of the foregoing results a shell with middle surface equation 


z= Cy + 3e,2" + 3e.y’ (22) 
simply supported along edges x = +a and y = +b. We further assume that the sup- 
porting structures are rigid in their own planes but non-resistant perpendicular to their 


own planes. 
The boundary conditions for w are 


Ow 


x= +0; w= 0, 3 = 0, (23a) 
Ox 
a . 
y = +b; w= 0, —s = (). (23b) 
dy 
The boundary conditions for F can be shown to be 
oF oF 
x= +4; : 4 = (), ——- = (0), (24a) 
Oy Ox 
= Al a 7 
y = +b; oF = Q, oF = (0. (24b) 
. Ox Oy 


The differential equations (I) and (II) reduce to the following form 


V’V'F = —Bite = + ¢ eu), (25a) 
oy Ox 
roy 2, OF, OF ™ 
DV Vw = —ph aE +c ay? + ct ax (25b) 


The form of these equations and of the boundary conditions (23) and (24) allows 
the following modes of vibration. 


) cos Oons1% COS BonsiY 


(26a) 


w= Ae ml 
SIN a@,X SIN Bo,Y 
\ COS Qon+1% COS BonsiY 
F = Be'**< (26b) 


sin Qo,t SIN Bo,Y 
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where 


mr mar - 
On = >) » = >,° (27) 
2a 2b 


Introduction of (26) into (25) gives 
Eh(c,B, + c20,)A,, — (a, + 82)’B,, = 0, (28a) 
[Diaz + B32)” — phwrJAr, + (c:8% + c20%)B,, = 0, (28b) 


where p and g may be any integers greater than or equal to unity. 
Vanishing of the determinant of (28) leads to a frequency equation of the form 


— oT n2 , EB (exo, + ¢,87)\" 
We = — (a, + 8 + —(|-—F-— !. (29) 
ph p\ Ho t+B 
Spherical cap with quadratic base. For a more detailed discussion of the frequency 
equation (29) we choose a shell for which 
—H —H 


a a 


a= b. (30) 


o=H, «a= 


Introduction of (30) into (29) transforms (29) into the following form 


4/2 2\2 , 2 9?) ao 2 2 
tT(p +q) Eh E i 192(1 — vr) HI. (31) 


192(1 — »°) pa* rip+q) h’ 


Woe = 


The frequency w,, is smallest when p = gq = 1. The value of w,, is given by, 
r E h’ 48(1 — v°) H” — 
we i =F) pa E +r T a Ss?) 
It is seen that the shell-curvature correction term is relatively largest for the lowest 
mode and becomes significant as soon as H/h is about unity. 

The frequency equation (31) may be compared with the relations (18a) which were 
involved in the justification of the use of (I) and (II). We have here both ZL and 1 of 
order of magnitude a. Accordingly, w as given by (31) is included in the range of admis- 
sible values which is determined by (18a, i) and (18a, ii). 

The following further observation may be made. As H/h increases the order of mag- 
nitude of the frequencies of transverse vibrations increases from (E/p)'’*(h/a’) to 
(E/p)'”(H/a’). On the other hand the order of magnitude of the frequencies of longi- 
tudinal vibrations of flat plates is (E/p)'’*(1/a). Since H/a < 1 for shallow shells, we 
see that, while shell curvature increases frequencies of transverse vibrations in relation 
to the corresponding flat plate frequencies, these shell frequencies remain of a lower 
order than longitudinal flat-plate vibration frequencies. It would seem that this fact is 
the physical reason why the transverse vibrations of shallow shells take place with 
negligible amount of longitudinal inertia. 
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TENSORS ASSOCIATED WITH TIME-DEPENDENT STRESS* 


BY 
BARBARA A. COTTER AND R. S. RIVLIN 


Brown University 


Abstract. It is assumed that six functional relations exist between the components of 
stress and their first m material time derivatives and the gradients of displacement, 
velocity, acceleration, second acceleration, --- , (n — 1)th acceleration. It is shown 
that these relations may then be expressed as relations between the components of 
m + n + 2 symmetric tensors if n > m, and 2m + 2 symmetric tensors if m > n. Ex- 
pressions for these tensors are obtained. 

1. Introduction. It has been shown by Rivlin and Ericksen [1]? that if we assume 
that the components of stress ¢t;; , in a rectangular Cartesian coordinate system 27; , 
at any point of a body of isotropic material undergoing deformation, are single-valued 
functions of the gradients of displacement, velocity, acceleration, --- , (n — 1)th accel- 
eration in the coordinate system x; at the point of the body considered, then the stress 
components ¢;; may be expressed as functions of the components of (n + 1) symmetric 
tensors defined in terms of these gradients. 

We describe the deformation by 


x, = 2(X; , b, (1.1) 


i 


where x; denotes the position, at time ¢, in the coordinate system x; , of a material 
particle which was located at X,; in the same coordinate system at some other instant 


of time 7. Let »;"’, vi, v{°, --+ , v;" denote the components of velocity, acceleration, 
second acceleration, --- , (n — 1)th acceleration at time ¢, in the coordinate system 7, , 
of a material particle located at x; . Then, if we assume 
dz, dv," dv,” av,” 
tas = tH, —, ee ee ee] (1.2) 
' aX,’ dr,’ dx, ax, 
it follows [1, Sec. 15] that ¢;; may be expressed as single-valued functions of the com- 
ponents C;; , A,; (r = 1, 2, --- , n) of (mn + 1) tensors defined by 
C — Ox, Ox; bey pas dv," dv;"’ 
“*  @X, dX,’ sith On; Ox; 
and 
DA;*’ av,” avs” 
a ae ind? P hid —m 4 a: a 1.3 
’ Dt sa " Oz; Ox; (1.3) 


where D/Dt denotes the material time derivative. This result was obtained from the 
consideration that the form of the dependence of the stress components ¢;; on the gradi- 
ents of the displacement, velocity, acceleration, --- , (n — 1)th acceleration must be 
independent of the particular choice of the rectangular Cartesian coordinate system z; . 

It will be shown in the present paper, from similar considerations, that if, instead 


*Received Aug. 6, 1954. The results presented in this paper were obtained in the course of research 
conducted under Contract N7onr-35801 between Brown University and the Office of Naval Research. 
tNumbers in square brackets refer to bibliography at the end of the paper. 
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of describing the dependence of the stress components on the deformation by six func- 
tional relations of the type (1.2), we have six independent functional relations of the 


form 
Iz, Ov, dv,” dv,” DM Pt D"t,. 
f,( 2% ; > ’ < —_, — = a pa ) at og “7 r; = ee ‘ne) = 0, (1.4) 
OX,’ Or, O2, Ox, pt” Di Dt 
with 
Fi _ fii 5] (1.5) 
then these functional relations must be expressible in the form 
Pa) A ee, + 2 4 Ee ess ot (1.6) 
if n > m, and in the form 
Patash fee se 4 ee eo Ee ee: (1.7) 
if n < m, where F;; = F;, in both cases, C,, and A,;’(r = 1, 2, --- , n) are defined by 
(1.3) and B,;’(r = 1, 2, --- , m) are the components of symmetric tensors, defined by 
ie r-1) OV; dv; 
BY? = —— + BY? Bo? — 
vs Dt dinate OX; 7 OX; 
and 
Bi; = t;. (1.8) 


Zaremba [2] introduced a rate of change of stress tensor, which is given in terms of 


the tensors ¢;; , A‘;’ and B};’ by Bi; — 4t,;Ai;’ — 3t:,Aj;. 


It may be remarked that Eqs. (1.4) and (1.5) are not, in general, sufficient for the 
determination of the stress resulting from the subjection of the material to a specified 
deformation history. They may be regarded as a set of six independent differential 
equations in six dependent variables ¢;;(¢;; = ¢;;) and one independent variable ¢. Suit- 
able “‘initial’”’ conditions at specified values of ¢ must be chosen if the equations are to 
have a solution. However, we are concerned here only with the limitations which must 
exist on the form of the relations (1.4), as a result of the necessity that they are invariant 
under a transformation from one orthogonal coordinate system to another, quite apart 
from any question of the sufficiency of the equations for the determination of the stress 
components. 

2. The deformation tensors. It is well known that if ds is the distance at time ¢ 
between two material particles of a body, undergoing a deformation described by (1.1), 
which are located at x; and x; + dz; in the rectangular Cartesian coordinate system 7, , 
then 

(ds)? = dz, dx, (24) 


Ox, OX; 

OX; 0X; 
where X; and X; + dX; are the positions of the particles at a previous instant of time 7’. 
Differentiating (ds)* r times with respect to t, we have, from (2.1), 


D'(ds)" ge le 
= A;; dx; dx; , ( 


dX; dX; , (2.2) 


bo 
(eS) 





bo 
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where A‘; is a symmetric tensor given by (1.3). A corresponding result in a convected 
coordinate system was obtained by Oldroyd [3]. 

Equations (2.2) and (2.3), with the left-hand sides given constant values, describe 
the deformation quadrics at the point of the body considered. 

It has been seen [1, Sec. 10] that A;;’ may also be expressed as 


av," + re x > ( r) ” dvi” 
Af? =* 2. 
Ox; + p/ Ox; a2; GA 
3. The stress tensors. If we define a sianlie ® by 
@ = t,; dz dz; , 3.1) 


then, since ¢,;; transforms as a tensor from the rectangular Cartesian coordinate system 
x, to any other, irrespective of any relative motion of the two coordinate systems, it is 
apparent that ® is a scalar, invariant under such transformations of the reference system. 
From (3.1), we have 


D'@ D*(t;; dx; dx;) Dt, D*(dz; dx;) 
wD ted «3 (1) Doty DM 





Dt Dt’ 2\o) Di Dt 
_s Dt; ¥ (4 ees D ae | 9 
= >| (3); pe 2\i) pe de | 6.2) 
Since 
D'(dx;) a) avs” 
7 ule dv; = ax, dx, 
and 
D*"'(d: nin a 
we — Ri = dy?” oo dx, , (3.3) 
where v‘° = z, , so that dv; /dz; = 5;; , we obtain from (3.2), 
D'® — I} (s\ Dt; a) é dvi" dvs” 
De & (8) Dt x t a, oe, ata 


With the definition 


: : s\ D’~*t,, << (q\ avi*~” av,” 
B *) - (‘): 2h ( ) ve a 5 
” 2, | q Dt? do i} @s; 2d; J’ (3.5) 
we can re-write (3.4) as 

‘® 


De = BSD dx, dz; , (3.6) 


with BS? = B;?. Since D'@/Dt' transforms as a scalar, between two rectangular Car- 
tesian coordinate systems with arbitrary relative motion, B{;’ transforms as a tensor 
between such coordinate systems. 

We note, from (3.6) and (3.3), that 


a+l1 
Bi;*" dz; dz; = aan = Di D (pw dx; dx;) 
3.7 
DB‘” ov’? yi? ( ) 
= (PB 4 By & — + BY S ) ae, de, 
dt ax, 
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Whence, 
+1) _ DBii (s) dv," Stel dv; 
By 3 HE + Re + ae (3.8) 


4. The stress-deformation relations. We assume that the dependence of the stress 
components on the deformation is described by the six functional relations (1.4) and 
the forms of the functions f;; are independent of the rectangular Cartesian coordinate 
system in which Eqs. (1.4) are expressed. Let x* be a rectangular Cartesian coordinate 
system moving in an arbitrary manner with respect to z,; and related to x; by 


2% = a;,;(x; — b;) a;;.% = 6%, (4.1) 


where a;; and b,; are, in general, functions of time. Let X* denote the coordinates in 
the system x* of a point located at X;, in the coordinate system x, and let v¥’, v¥, 
--- . v*" be the components of the velocity, acceleration, --- , (n — 1)th acceleration 


respectively in the coordinate system x*¥ . Then, if ¢* are the components of the stress 


in the coordinate system x* , we have 


ac* avt? av*? ov*'” Dt* Dt ‘) 
f (: : p Soo ae Ge SM... =) oe O, (4.2 
: re Dy Dt” ' 4.2) 


OL 


~A5 ea 


It has been shown in a previous paper [1, Secs. 5 and 15] that we can choose the 


ox* 


coordinate system x* in such a way that: 


(i) instantaneously at time ¢, the directions of the axes of x* are parallel to those 
“ ? 

of x; , so that 

a;; = 6; ; (4.3) 
the angular velocity, angular acceleration, angular 
second acceleration, --- , angular (n — 1)th acceleration of the coordinate 
system x* relative to x; are such that the velocity, acceleration, , (n — 1)th 
to the coordinate system x* , in the immediate 


(ii) instantaneously at time ¢, 


acceleration fields relative 
neighborhood of the material particle considered, are irrotational; 
instantaneously at time 7’, the axes of the coordinate system x* have directions 


(iii) 
relative to the coordinate system x; defined by 

OX, = 

Gi; =ay (OC de, (4.4) 
i ax k 
where 

Ox; OX v 
C = || sy-ar (4.5) 

|| OX, OX, 


and c(= || ¢,; ||) is the matrix satisfying this equation which has all its eigen- 
values positive. 
The choice of the coordinate system x* in accordance with the condition (ii) implies 


that, at time ¢, 
(4.6) 


ov**"’ /dz¥ = dv*""’ /dx* (ry = 1,2, --- ,n). 


Also, the choice of the coordinate system x% in accordance with conditions (i) and (iii) 
implies that, at time ?, 
62% /OXF = ¢;; . (4.7) 





to 
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With the notation 


\ 


oy — 1 fav , avg 
dy = 2 art + “On* ), (4.8) 





it follows from (4.6) and (4.7) that Eq. (4.2) can be re-written as 


Dig... Dia) _ 4 
Dt’ a dt 








(4.9) 


* * (1) (2) — * (n) * 
J Se, ; rl ) = is ’ dy, ’ tye ’ 


at time ¢. 
From (3.5), (4.6), and (4.8), bearing in mind that v*’ = x* and t* = t* , we see that 
if B*“? are the components of the tensor BS’ in the coordinate system 2x* 
i p i o ’ 


oF D'-t3 < (2) av’ ot] 
* - _foh ak 
Bx? = >| (’) Di 2 Ai) axt ort 


q@=0 


¥ ud qa a 
~ Beh 1 (9) Ses E (Qareran” 
“2 Y 1=0 


@ 





(4.10) 


Also, at time ¢, from (2.4), (4.6) and (4.8), we see that if A*;*’ are the components of 
the tensor A;;’ in the coordinate system x* , then at time ¢, 


At," = 2d¥," + Oi "Yaseen (4.11) 
Since, 

At = Ode, (4.12) 
we see, from (4.11), that d*"’ can be expressed as a polynomial in the quantities A*” 
A*® ... , A**”, We also see, from (4.10), that D*t*/Dt' can be expressed as a polynomial 
in the quantities A*,”’, A*,”, --- , AX”, t% , BR”, BR”, --- , B*?. Consequently, 

Eqs. (4.9) may be re-written in the form 
Cle» A”, AR <** AR, aR AR ++ (Fo G, (4.13) 


if nm > mand in the form 


P:iCre » A, AN, -°* A”, te Be, BR, -** BR”) = 0, 4.10 


if m > n. It may be noted that if, in (4.9), f,;; is a polynomial function of c,, , d%{"’ , 
, D"t*/ Dt", then in (4.13) and (4.14), ¢,;; is a polynomial in the dependent variables. 
Since B;;’ and B*;"’ are the components of the same tensor in the coordinate systems 


x, and x* respectively, we have 


Bor’ = Bi; Gpile and t% = t;;Q5iBq; « (4.15) 


Since A};’ and A*;”’ are the components of the same tensor in the coordinate systems 
x, and x* respectively, we have 


ya 


Pa 


) me As; Gye; - (4.16) 


Since the coordinate system is chosen in accordance with condition (i) we see that, 
at the instant of time ¢, a;; is given by (4.3) and Eqs. (4.15) and (4.16) yield 


Be? = B® and At”? = AM, (4.17) 








182 BARBARA A. COTTER AND R. S. RIVLIN [Vol. XIII, No. 2 


~ 


Introducing the results (4.17) into Eqs. (4.13) and (4.14), we have 


weit | ee, ae, **? oes er ee ee, *** 5 ES oe Be (4.18) 
ifn > m and 
ee = eo dor 4 *** . oe > es ge ee oc *** oe ee Oe (4.19) 


if m > n. 
Following the method adopted by Rivlin and Ericksen [1, Sees. 7 and 15] and em- 
ploying the notation 


Ox; O2, 


al ee a, on 
the relations (4.18) and (4.19) may be written as 
Pie eee ae os, Bk Be ss oe @ (4.21) 
if n > m and 
EE bo Ee bra Se i ay OS 1.22) 


if m > n, where F,; is a single-valued function of the independent variables. 

If we assume in Eq. (1.4) that the functions f;; are polynomials in the variables, 
then it follows from the manner in which Eqs. (4.21) and (4.22) are derived that F;,; are 
polynomials in the variables. It is also readily seen that if we assume the functions f,; 


are single-valued functions of dv,"’/dx, , dv,” /dx, , +--+ , dv," /dx, , ty, , Dt,,/Dt, D°t,,/Dt’, 


-++,D"t,,/Dt” then they may be expressed in the form (4.21) or (4.22) with C,, omitted. 
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STRESS ANALYSIS IN VISCO-ELASTIC BODIES* 


BY 
EK. H. LEE 


Brown University 


Abstract. The analysis of stress and strain in linear visco-elastic bodies is considered 
when the loading is quasi-static so that inertia forces due to the deformation are negli- 
gible. It is shown that removal of the time variable by applying the Laplace transform 
enables the solution to be obtained in terms of an associated elastic problem. Thus the 
extensive literature in the theory of elasticity can be utilized in visco-elastic analysis. 
The operation of the transform on the prescribed boundary tractions and displacements 
and body forces may completely modify the spatial distribution in the associated 
problem. For proportional loading, in which the space and time variations of the pre- 
scribed quantities separate, the spatial distribution is maintained in the associated 
problem. A convenient method of treating a common case of non-proportional loading, 
moving surface tractions, is demonstrated. This work is compared with related approaches 
to this problem in the literature of visco-elastic stress analysis. 

1. Introduction. We are concerned with quasi-static problems in which inertia 
forces due to the deformation are negligible. Inertia forces due to effectively rigid body 
motion such as centrifugal forces come within the scope of the analysis. We shall consider 
materials satisfying the general isotropic linear visco-elastic law: 


Ps;; ins Qe; ; ? (la) 
P'a;; = Q'e:: , (1b) 


where P, Q, P’ and Q’ are linear operators of the form >°% a,D", and D is the time de- 
rivative 0/dt. The coefficients a, and the numbers m are in general different for each 
operator, although certain restrictions on the m’s are required to determine observed 
physical characteristics. ¢;; and e,;; are the stress and strain tensors, the latter considered 
to be infinitesimal, and s;; and e;; are respectively their deviators defined in the usual 
way: 


. — 1 sil —s 
me Rie 37 Kk OG; ’ Cis = i; 3€kk Oi; ’ (2) 


where 6,; is the Kronecker delta 6;; = 1,7 = 7; 6,; = 0,7 #/. 
Equations (1) are a generalisation for combined stresses of the well-known visco- 


elastic law, 
Po = Qe, (3) 


relating stress and strain when only one component of each tensor is needed, as for 
example in simple tension (see Alfrey [1]}). In developing the most general linear isotropic 
relation, operator equations of the type (3) can be written between the deviators and 
the first invariants of o,; and e,;; as detailed in (1). This is analogous to the expression of 
the general isotropic elastic relation in terms of the two constants, shear modulus and 
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bulk modulus. (1) is equivalent to the relations used by Read [2], and is more general 
than that prescribed by Tsien [3] according to the same premises. 

We adopt the form (1), which corresponds to a discrete spectrum of relaxation times, 
in preference to the hereditary integral method of specifying visco-elastic behaviour, 
which gives a continuous spectrum, because of simplicity and since we are particularly 
interested in stress analysis for short loading times. In such a restricted range visco- 
elastic behaviour can be represented with tolerable accuracy by relations of the form 
(1) with only a few terms in each operator. This is equivalent to representation by : 
simple model of dashpots and springs. Methods of determining such simple models to 
cover a limited frequency range will be described elsewhere. The hereditary integral 
representation equivalent to (1) has been used by Volterra [4]. 

2. The equivalent elastic problem. As shown in Fig. 1 we consider bodies subjected 


Fia. 1 


to prescribed body forces f,;(x; , 4) per unit volume, tractions 7,(z; , t) on the surface S 
or surface displacements u,(z; , t), including combinations of these in the groupings in 
which they occur also in elasticity theory. In order to determine a complete solution we 
must obtain the stresses ¢,;(2, , t) and the displacements u;(z; , t) throughout the body 
V to satisfy the stress-strain relations (1), and the equilibrium equations: 


6:54 = Su, d, (4) 


where as usual the subscript after the comma indicates differentiation with respect to 
the corresponding space coordinate. The strain tensor is given in terms of the displace- 
ment by 


é; = d(u jo U;,;). (5) 


The displacement must be compatible with the prescribed surface displacement, 
and the stresses with the prescribed surface tractions according to: 


T; = o,;n; on S, (6) 


where n, is the outward normal. 

We now remove the time dependence by operating with the Laplace transform on 
all these equations (see for example Carslaw and Jaeger [5]). The transform is denoted 
by a star on the corresponding function. We will consider problems in which the bodies 
are initially undisturbed since this is the common situation, but the Laplace transform 
technique is also well suited to problems with non-zero initial conditions. With zero 
initial conditions an operator in D merely becomes the same function of the transform 
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parameter p, so that the required equations become 


P(p)st; = Q(p)et; , (7a) 
P'(p)ot; = Q'(pe% , (7b) 
o*, ; = ft(u ,p), (8) 
ef, = 3(u?; + u?.,), (9) 
Tx, ,p) = etm; , (10) 


(7)-(10) represents a stress analysis problem for an elastic body of the same shape as the 
visco-elastic body with elastic constants a function of the parameter p. The prescribed 
surface tractions T*(x;, p), body force f*(x,; , p) and boundary displacement u*(zx; , p) 
are also functions of p. When the stresses o*(a, , p) have been determined, inversion of 
the Laplace transform gives the desired stresses o;;(x, , t) for the visco-elastic problem. 
and similarly for other variables. Thus we can utilize the extensive literature of the 
theory of elasticity to evaluate visco-elastic stress analysis problems. 

Although the body shape for the associated elastic problem is the same, the prescribed 
distribution of surface tractions, displacements, and body forces may be quite different 
since the Laplace transform ¢*(x,; , p) will in general have an entirely different space 
distribution than g(x; , t). Thus the associated elastic problem is not in general simply 
related to the visco-elastic problem. There is however a common type of problem, which 
we shall call proportional loading, in which the space and time dependence of prescribed 
quantities separates out with a common time dependence. Thus 


T(x, ,t) = Tix) fd, 
fix; ,0 = fiz) f(d, (11) 
ul(a,;) f(b. 


u(x; , t) 


In this case the Laplace transform simply changes f(t) to f*(p) leaving the space 
dependence unchanged. Since f*(p) then merely appears as a multiplying factor, the 
associated elastic problem contains the same spatial distribution of prescribed quantities 
as the visco-elastic problem. This represents a considerable simplification since the 
associated elastic problem can be solved for T/ , f; and u{ prescribed, and these are 
independent of p. The analysis being linear, if the loading falls into groups of proportional 
loading and prescribed displacements each with different time dependence, they can be 
considered separately and superposed. An important class of problems which does not 
fall in this category is associated with moving loads. A particular method of treating 
these is detailed below. 

3. A simple example. Let us consider the problem shown in Fig. 2 of a vertical 
point force P(t) acting normally at a fixed point on the surface of a semi-infinite visco- 
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elastic body. The point of action is the origin of cylindrical coordinates (r, 0, z), with 
the body occupying the space z > 0. Since only one point force is specified this is a case 
of proportional loading, and the associated elastic problem is that of a semi-infinite 
elastic body with a normal point load of magnitude P*(p). The solution to this elastic 
problem is given by Timoshenko ([6], p. 331), and to demonstrate the method we will 
consider the stress component o* which is given by 


of = — 4{(1 — 2»)| 4 _ 5 (r? + 2°)7" ‘| — 3r2r° + zy (12) 


where vy is the Poisson’s ratio of the elastic body. Transforming the notation of (7) to 
more familiar elastic constants 


Q(p)/P(p) = 2G, (13) 
where G is the shear modulus, and 

Q’'(p)/P'(p) = 3K, (14) 
where K is the bulk modulus, since o* is three times the average hydrostatic tension. 
We can express Poisson’s ratio in terms of K and G (see [6] p. 10) according to the 
relation 


3G 
= ra F 
ides es said 
Thus, in the notation of (7), (12) becomes 

P*(p) 3Q(p)/P(p) E € us — | se 5/2 

t= SSC 42) | - Bre? + 2)°”?. 

" 2r \[2Q0'(p)/P (p)] + [(Q(p)/P@)| Lr r - _ J 
(16) 


For a prescribed material, with known operators P, Q, P’ and Q’, the inverse of (16) 
gives the radial stress in the visco-elastic problem. Suppose, for example, that in shear 
the material behaves as a delayed elastic Voigt material, but is perfectly elastic under 
hydrostatic pressure. (1) would then read 


(AD + Boe;, , (17a) 


va) 
I 


oi: = Ce; , (17b) 


where A, B and C are constants of the material. Suppose that a constant load P, is 
suddenly applied at ¢ = 0, and maintained. Then P(t) = P, ,t > 0, and P*(p) = P,/p, the 
transform of the Heaviside step function (see [5] p. 4). Thus (16) becomes: 


P, f 3(Ap + B) ] ee —— a Pye a 
” = Dep \2C + Ap +B Pp 2? | — Bre + 2)). (18 
C; 2rp \2C + Ap + ml a y+) | ar°z(r° + 2°) J (18) 


r 


This can be inverted directly by partial fractions to give: 
oe = — \aq- a ([B UH C2C exp [—(2C + Be Ail 4 =s (r° +2) 1/ | 


— 3re(r? + ie t> 0. (19) 
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Thus for the visco-elastic problem with a constant force applied, the stress consists 
of two terms, one remaining constant, and the other dying out exponentially with time. 
There is thus a relaxation influence between the shear and the hydrostatic response of 
the material. 

For more complicated materials, Q/P and Q’/P’ will still be rational functions of 
(p), so that this method of inversion can be used, although the manipulation will be more 
cumbersome. 

For more complicated P(t), it might be worthwhile to separate the inversion of 
P*(p) and of the function of p arising from the material properties, since these appear as a 
product. The inversion can then be expressed as a convolution integral (see [5] p. 7) for 
general P(t). Thus for the material considered in the examples, (18) can be inverted in 
this way, with the special value P,/p replaced by P(p), to give: 


1 fr 2C ' ; 
to = | 3 0,7) - 7 exp [—(2C + B)x/A\H() [PC — dr 
(20) 
x E “Ss (r? + 27)" ‘| — P(t)3re2(r? + 2°)? : 
; ; 
where 4(0, 7) is the Dirac delta function, and H(r) the Heaviside step function. For the 
special case P(t) = P, , the integral is simply evaluated to give (19). 


4. Moving loads. A common type of loading which falls outside the category of 
proportional loading is that due to moving surface tractions. Such a situation is often 
produced mechanically or by the pressure system in a moving fluid. It can be handled 
directly as described in the previous section by taking the Laplace transform of the 
moving surface tractions, but an alternative method is described which may be more 
convenient. 

Instead of taking the transform of the prescribed boundary conditions and body 
forces to determine the associated elastic problem, and then solving this, the prescribed 
visco-elastic problem could first be solved assuming it to be elastic, and the transform 
of the resulting stresses would be the desired solution of the associated elastic problem. 
However the surface tractions vary in time, at each instant the solution of a quasi-static 
elastic problem is simply the solution of the problem with the current tractions held 
constant. This is of course not true for a visco-elastic problem in which the history of 
loading has an influence, as was demonstrated by the example of the previous section. 
Thus, the solution of the elastic problem for the series of instantaneous distributions of 
surface tractions gives the sequence of stress values at each point, and the transform of 
these values gives the required solution of the associated elastic problem. The elastic 
constants are the functions of the transform parameter p associated with the trans- 
formed stress-strain relation (7), but since they are not time dependent this adds no 
complication to the determination of the transform of the stresses. On inverting to 
obtain the stress for the visco-elastic problem, the fact that these elastic constants are 
functions of p modifies the reverse process. 

To illustrate the method, let us consider the problem, shown in Fig. 3, of a point 
load P(t) moving along the x axis on the surface of a semi-infinite visco-elastic body. 
The motion of the force is prescribed by the displacement from the origin &(¢), and the 
force is first applied at — = 0, ¢ = 0 with the body initially undisturbed. At each instant 
the solution of the elastic problem with the same loading is that used in the previous 
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section, it merely being necessary to change the origin of coordinates. To obtain a simple 
example for illustrative purposes, let us consider the stress component c, at the point 
(xz, 0, z) directly under the path of the load. At any instant the stress of in the elastic 


problem with this point force loading is given by (12), with r replaced by [x — £(¢)], 

ee ee 
ce: 
5 aes yy 

+< ———— a i 

\ f / 

\ ; ff 

Fia. 3 


since the directions of x and r are the same at the point in question, and the radius of 


the point (x, 0, z) relative to cylindrical coordinates based on the load axis is [x — £(t)]. 
Thus: 
. Pas 1 2 , oe 
a A — 2p) a - —<3 (a — &) +2} - 
2r (2 — &) (x — &) 
21) 


— Be — sale — BF + 2 


where » is the Poisson’s ratio of the elastic material. o*(zx, z, p) of the associated elastic 
solution is the Laplace transform of this, with 1 — 2» replaced by the transform of the 
equivalent visco-elastic operator [3Q(p)/P(p)]/[2Q’(p)/P’(p) + Q(p)/P(p)], as in 
Eq. (16). In carrying out the transform the ¢ appearing in the prescribed £(¢) must be 
included. The stress c,(z, z, t) in the original visco-elastic problem is now obtained by 
taking the inverse of c% which contains the transform parameter p both from the trans- 
form of c: (21) and from the visco-elastic operators. The transform of (21) may be 
quite complicated depending on the particular forms of P(t) and &(t). However since the 
elastic constant v appears only as a multiplying factor there is no need to carry out the 
transform of oc , since the inverse transform can be effected using the product rule in 
terms of a convolution integral (see [5], p. 7). This gives for the stress component a, at 
the point (z, 0, z) for the particular visco-elastic material considered in the previous 
section the value: 


1 {ft 2C 
o,(z,2, 0) = on | | 3[6(0, 7) — 4 OXP [—(2C + B)r/A]H(7)|P(t — 7) 
l Z if aii 2 2)-1/2 (99 
P= ol” G@eo ar On eee fe aaa 
le \ 7 . $ q 


— P(H3[x — e(HPeh[e — EOP + 24°? 


In many cases it may be much easier to complete the solution in this form, in place 
of carrying out the Laplace transform of 7';(x; , ¢) to obtain o*(x,p). In using the method 
of the present section, the need for determining this quantity is avoided. 

Although this alternative method is particularly well suited to a problem with 
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moving tractions which apart from this motion are otherwise unchanged in form, it can 
also be used for the general case 7';(z; , ¢). In this case a series of elastic solutions would 
be needed, with the space distribution of traction that obtaining in the visco-elastic 
problem at each instant. 

5. Discussion. The analysis given in this paper was developed as an extension of 
the work of Alfrey [7] and Tsien [3]. Alfrey showed that for an incompressible visco-elastic 
medium with prescribed surface tractions the stress is the same as that for an elastic 
body. This result is obtainable from the present analysis by taking P’/Q’ equal to zero. 
Tsien showed the same to be true for a compressible visco-elastic body under proportional 
loading, but he had assumed a restrictive relation between the operators P, Q, P’ and Q’ 
such that the equivalent Poisson’s ratio reduced to a constant rather than a rate de- 
pendent operator. The example of Sec. 3 shows that in the general case the stress field 
changes under constant tractions, so that it cannot be equal to the solution of an invariant 
elastic problem. 

In the work of Alfrey and Tsien only a single pair of operators P, Q appeared, and 
it was possible simply to separate the time and space operators in the differential equa- 
tions, and to obtain the elastic solution from the space operators. For the general isotropic 
medium four rate operators P, Q, P’, Q’ arise, and the Laplace transform assists in 
separating out the time dependence. The first application of a method of this type 
which I have seen was the use of Heaviside’s operational calculus for the stress analysis 
of certain visco-elastic bodies by Jeffreys [8]. More recently Read [2] has used the Fourier 
transform to discuss the general dynamic problem. He obtains an associated elastic 
problem similar to that developed in Sec. 2 of the present report, but from the standpoint 
of application this result is extremely restricted in the dynamic case, since the inertia 
forces lead to body forces in the associated elastic problem which are proportional to 
displacements. Solutions of this type of problem are not common in the elasticity litera- 
ture. For the quasi-static case the principle of the present paper is closely related to 
Read’s, but his illustration of the use of the operator method is restricted to the case of 
proportional loading. For the usual type of problem in which loads are applied with or 
without an initial disturbance in the medium, it is felt that the Laplace transform is 
more convenient to use than the Fourier transform. The examples given in this report 
indicate the convenient form in which solutions are derived. 

Papers by Mindlin [9] and Graffi [10, 11] are concerned with the special situation 
when the visco-elastic stress field is identical with an elastic field. The former studies 
this question from the standpoint of photo-elastic testing. Special solutions in which the 
space and time variables separate have been considered by Volterra [4]. The example of 
Sec. 3 shows that this condition does not obtain for a constant fixed force, and so is 
quite limited for applications. Another particular group of problems which has received 
considerable attention is when the time dependence is entirely represented by a factor 
e‘**. When this is factored out, an equation in space coordinates only is obtained. The 
relation of this type of problem to both elastic and compressible-viscous fluid theory is 
stated by Oestreicher [12], and an example is given. 

This assessment of the literature suggests that the generality of the analysis pre- 
sented in this paper is needed for many problems of visco-elastic stress analysis, and 
that the use of the Laplace transform is a convenient means of treating such problems. 
The ease with which solutions can be evaluated depends on the way in which the elastic 
constants appear in the solution of the associated elastic problem. Much is known about 
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this, for example, as stated by Read [2], that they do not appear in the expressions for 
stress in the plane problem with prescribed surface tractions, if the tractions on internal 
boundaries are in equilibrium. In most elementary solutions they appear in a simple 
form so that the use of these solutions in visco-elastic analysis leads to comparatively 


simple evaluations. 
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—NOTES— 


NOTE ON A NON-HOLONOMIC SYSTEM* 
By O. BOTTEMA (Technische Hogeschool, Delft, Netherlands) 


A mechanical system with the coordinates g, , g2 , «++ qd, is called non-holonomic if 
there exist one or more non-integrable kinematical relations of the type 


A, dq, + A, dq2 + re + A, dqn _ 0, (1) 
where A; is a function of g; , G2 , --+ , g, - If the number of relations is m the system 
has n — m degrees of freedom. 


In every textbook on mechanics which deals with the matter examples of non- 
holonomic systems are given. We always meet the sphere moving on a rough plane: 
the conditions which express that the velocity of the point of contact is zero are non- 
holonomic (n = 5, m = 2). A more simple and very attractive example has been given 
by Carathéodory’ in a paper on the sleigh: a line-segment AB is moving on a horizontal 
plane and one of its points C is subjected to the condition that its velocity has no com- 
ponents perpendicular to AB(n = 3, m = 1). In all these cases the non-holonomic 
constraint is due to friction. 

An example of a non-holonomic system of another type seems the following. Consider 
a (horizontal) dise D, which is able to rotate without friction about a vertical axis / 
intersecting D in O, and a particle P which moves without friction on the disc. If there 
are no external forces acting on the system its moment of momentum about / is constant. 
We add the condition that this constant is zero; the system is now seen to be a non-holonomic 


one (n = 3, m = 1). 


If OA is a fixed line on the disc, OX a fixed horizontal line in space, XOA = g, 
OP =r, AOP = yy, then r, ¢g and y are coordinates of the system. J being the moment of 
inertia of D about the axis /, and m the mass of the particle P, the condition reads 


(I + mr’) de + mr dy = 0. (2) 


Obviously this relation is non-integrable. It can also be shown directly that there 
cannot exist a relation between the three coordinates. Suppose that according to a force 
acting between the particle and the disc, P describes a curve r(t), Y(t) on D, starting 
from rest in A(r, , ¥,) and moving to B(r, , ¥2). Then during this relative motion the 
dise rotates about an angle 

ate 2 
g=-m | Fi rr y(t) dt. 
If B coincides with A this value is in general not zero: that means that the coordinate ¢ 
is not determined by r and y. We have a very simple case if mr* can be neglected with 
respect to J; then ¢ = —m f;? r°y’ dt = —2mF, where F is the area enclosed by the 
curve which P has described on the disc. It is therefore possible to move P in such a 
way that to given values of r and y an arbitrary value of ¢ belongs. Suppose that the 
force acting between D and P depends only on r and y and that moreover this field 


*Received July 23, 1954. 
\Carathéodory, Der Schlitten, Z. angew. Math. Mech. 13, 71-76 (1933). 
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a“ 


force is conservative, the potential function being mV(r, y). Then the Lagrangian 
function i 
L=T—V =I + m’)e? + mre + 4m’? + 4mr? — mV(r, YW). 


™ 


The Lagrange equations of the second kind for ¢, y and r are respectively 


7 


a > ° re , > 
7 (I + mr)eo + mrp} = XJ + mr), 
— ig scl re 10V 
tr (ta a 2re te 2r y + ry a ae ae vr, 
r oy 


- 2 ve 2 av 
T — Fo _— 2re y ame ry - a mg = 0. 


From the first of these equations it follows that the multiplier \ is 0. g and yg’ can 
be eliminated and we have two equations for the relative motion r, y of the particle P 
on the disc. They are rather complicated. Of course we can take 7’ + V = c instead of 
one of them. If we have the case mentioned above, where mr* may be neglected when 
compared with J, the equations are 


‘ ‘Lge ae Ee 2 OV sce i — 
g=--—try, f=" + = eG; ry tr —2V =e. 
I or 

This means that the relative motion is the same as it would be if the disc were fixed. 
We give a simple example: V = k°/2(r° — 2ar cos W), the force on P being an attracting 
force directed to the center C(r = a, Y = 0), and proportional to the distance PC. If 
x=rcosy—a,y=rsiny, wehavexr’ = —k*2z,y = —k’y,x = C, cos kt + C, sin kt, 
y = C, cos kt + C, sin kt, where C; are constants of integration. Therefore 


r = (Cy - & cos’ kt + (cy + C3) sin kt — (CC, a C3C,) sin 2k¢ 
+ 2a(C, cos kt + C. sin kt) + a’, 


y = arctan 2m st a E wn 
C, coskt+ C,sinkt+a 


ry = ay + (ay — yx) = ay + k(C,C, — C.C)). 
Hence 
1a = k - \ Y 
on -—> (C, cos kt + C,sin kt) — - Ct —~ re 4+ €,. 


ON LINEAR INSTABILITY* 
By AUREL WINTNER (The Johns Hopkins University) 
1. Let the coefficient function of the linear differential equation 
x’ + fijx =0 (1) 


be real-valued and continuous for large positive ¢. Consider only those solutions x(t) 
of (1) which are real-valued and distinct from the trivial solution (= 0). Then, since 


*Received August 12, 1954. 
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the values of x and x’ at a given ¢ determine a solution of (1) uniquely, x and x’ cannot 
vanish simultaneously, and so z(t) must change sign whenever it vanishes. 

If there exists a f = 7’ beyond which the function z(t) does not change sign, so that 
a(t) * O for every ¢ > T, then (1) is called non-oscillatory, and otherwise oscillatory. 
This classification of the differential equations (1) is independent of the choice of the 
particular solution x(¢), since, in view of Sturm’s separation theorem, every solution of 
(1) must have zeros clustering at = © if one solution has such zeros. 

2. With reference to a differential equation of the form (1), three types of stability 
will be considered: (i) strict stability, requiring 

lim sup {2°(t) + 2”(t)} < © (2) 
for every solution and its derivative; (ii) bownded stability, requiring only the coordinate 


condition 
lim sup | z(t) | < @ (3) 


t-—@ 
for every solution; (iii) oscillatory stability, requiring 


lim N(i) = © (4) 
t-@ 
for some or for every solution, where N(¢) denotes the number of times the solution 
a(t) changes sign between a fixed ¢, and a variable t. 

That (i) is actually stricter than (ii) follows by choosing f(t) to be any function 
which tends to © sufficiently fast, with a sufficient degree of regularity in its growth, 
and possessing a continuous second derivative. For then it follows from a general asymp- 
totic formula’ that (3) holds for every solution but (2) does not hold for any solution. 

Dynamical considerations suggest that (iii) is necessary for (ii). Actually, not only 
(4) but even 

lim inf N(t)/t > 0 (5) 
must hold’ when (3) is satisfied by every solution. On the other hand, trivial examples 
show that (4) can be satisfied when (5) is not. Hence (iii) is necessary, but not sufficient, 
for (ii). 

If an arrow is meant in the sense of logical implication, then the preceding considera- 
tions can be summarized as follows: 

(i) — (ii) — (iil), (6) 
where neither of the arrows can be reversed. 

3. In what follows, f° will denote (as customary) the number max (0, f), that is, 
f =0if f Ss Oand f* =| f| if f = 0. Thus f* (4) is a continuous function, since f(é) is. 

The second of the implications (6) states that (1) must possess some solution satisfying 

lim sup | z(t) | = @ (7) 

t-—-@ 


if (1) is non-oscillatory. But the converse is not true (see the trivial example (9) below). 
It is therefore natural to raise the question for a condition which will allow (1) to be 


1Cf. A. Wintner, Phys. Rev. 72, 81-82 (1947). 
2A. Wintner, Quart. Appl. Math. 7, 115-117 (1949). 
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oscillatory but will bring (1) so ‘“‘close’’ to a non-oscillatory state that (7) will remain 
true for some solution. The character of such a condition on f(t) is suggested by the 
following consideration: 

Since (1) is non-oscillatory if f(t) = 0, it follows from Sturm’s comparison theorem 
that (1) is non-oscillatory if f(t) = 0, that is, if f*(t) = 0. Hence, what suggests itself 
is that (1) will have some solution satisfying (7) whenever f*(t) is “small” in some 


sense; say in the sense that 
either [ fi(jdt< @ or f(j—-0 (8) 


ast — o. It will turn out that either of the conditions (8), and much less, is sufficient 
to this end. 
Needless to say, if f*(#) is “very small,” then (1) can be non-oscillatory. But this 
will not be the case by virtue of (8). In fact, if 
f() = at”, (lst< 0), (9) 


where a is a real constant, then, while both conditions (8) are always satisfied, (1) will 


be non-oscillatory (if and) only if? 
—-~o <a s1/4(>0). 9 bis) 
4. Put 
F*(t) = f*(s) ds, (10) 
where /, is arbitrary and 4, S$ t < o. It will be shown that not only condition (8) but 


even the condition 
lim inf F*(t)/t = 0 (11) 
se 
is sufficient in order that (1) have some solution satisfying (7). 
In a certain sense, this criterion is the best possible of its type, since (11) cannot be 
relaxed to 
lim inf F’(t)/t = e, or even to lim F*(t)/t = e, 


t-0 t-2 


if « > 0. In fact, (10) shows that the latter condition is satisfied if f(t) is a positive 
constant (= e); in which case, however, (1) is stable in every sense. 

The proof of the italicized statement depends on two known facts. The first is the 
fact, quoted above, according to which the inequality (5) is a necessary condition for 
the bounded stability of (1). The second is the following fact*: If N(¢) and F*(é) are 
defined as above, then, without any assumption on (1), some constant multiple of the 


function 
1+ {iF*()}'? (12) 


will always exceed N(t), ast — . 
The latter estimate of N(¢) implies that 


lim sup N*(t)/{tF*())} < @. 


t9@ 


8This classical remark of A. Kneser follows by observing that substitution of z(t) = ¢® into the 


case (9) of (1) leads, for b = b(a), to a quadratic equation the roots of which are real if and only if 


fa S 1. 
4P, Hartman and A. Wintner, Amer. J. Math. 71, 206-213 (1949). 


1955] RICHARD BELLMAN 195 


Hence, if (11) is satisfied, then 
lim inf N*(é)/¢ = 0. 
t-@ 
Since this contradicts the necessary condition (5), the proof is complete. 
5. In view of the first of the conditions (8), it is worth mentioning what happens 
to each of the three stability properties, considered above, if (1) is replaced by a differ- 
ential equation 


2’ + gx =0 (13) 


which is a small perturbation of (1), in the following sense: 


[ 110 - Wo \a<e. (14) 


Under the assumption (14), both (1) and (13) are stable in the sense of definition 
(i) if either of them is, and both (1) and (13) are stable in the sense of definition (ii) if 
either of them is. Both of these criteria (which are quite independent, since (3) does not 
imply (2)) are contained, as corollaries, in known* asymptotic correspondences between 
the solutions of (1) and (13), when (14) is satisfied. The situation is changed if (i) or 
(ii) is replaced by (iii). In fact, if f(t) = t-*, then (1) is oscillatory (since (9) then holds 
for an a > 1/4), and (14) is satisfied if g(t) = 0, but (13) is then non-oscillatory. 


5A. Wintner, Amer. J. Math. 69, 261-265 (1947). 


ON PERTURBATION METHODS INVOLVING EXPANSIONS IN 
TERMS OF A PARAMETER* 


By RICHARD BELLMAN (Rand Corporation) 


Summary. It is shown by means of some examples from the theories of linear 
algebraic equations, linear integral equations and nonlinear differential equations that 
the effectiveness of the method of expanding a solution in a power series in terms of a 
parameter may in many cases be greatly increased by expanding in terms of a suitably 
chosen function of the parameter. This is particularly the case when the physical setting 
of the problem allows only positive values of the parameter to enter. 

1. Introduction. A standard tool in the theory of functional equations, of both linear 
and nonlinear character, is the expansion of the solution as a power series in a parameter 
appearing in the equation, or in the boundary conditions. Some typical examples of 
equations involving a parameter are 


(a) 2+2z=, 
N 
(b) a+r Lia wz=c, t=1,2,---,N, 
1 (1.1) 
(c) fa) +» [ K(e, IW) dy = g(a), 
(d) a’ + X(2” — 1)’ +2=0 
(e) v’ +2+dz* = 0. 


*Received Sept. 20, 1954. 
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Two possibilities arise when the formal expansion for the solution is obtained. The 
expansion may be an asymptotic series, convergent for no non-zero value of the param- 
eter, or it may be a power series possessing a non-zero radius of convergence. We shall 
be interested here only in the case where there is a finite radius of convergence, although 
the idea we present is equally applicable to the case of an infinite radius of convergence, 
in which case we may wish to speed up the convergence of the series, or to the case of 
an asymptotic series, in which case we may wish to improve the best possible approxi- 
mation. 

As we know, the radius of convergence of a power series is determined by the distance 
from the origin to the nearest singularity of the function. Consequently, even if we are 
concerned with determining the numerical value of a function in a region in which it 
has no singularities, we are very often prevented from using the power series for this 
purpose because of singularities which occur in regions of no interest to us. 

Let us give some examples which we shall discuss again below. If \ is small and positive 
the positive root of (1.la) has a power series expansion of the form 

1/2 
r=rA—-N+--- = oe tte. (1.2) 


_ 


This expansion may only be used up to the value \ = 1/4, because of a singularity at 


A= —-1/4. 
Similarly, if we consider the vector-matrix form of (1.1b) « + A\Az = ¢, the Liouville- 


Neumann solution 


x=c—rAc+::: (1.3) 
will have a finite radius of convergence determined by the location of the roots of the 
determinantal equation |/ + AA | = O. If we take A to be a positive definite matrix, 


and ) a positive number, there will always be a unique solution of (1.1b) for \ sufficiently 
large. Nevertheless, we cannot use the power series in (1.3) to determine the solution 
directly. 

The same remarks may be made concerning the linear integral equation in (1.1c). 

Turning to the differential equations above, the first of which is the famed equation 
of Van der Pol, and the second the normalized equation of a nonlinear spring, it is known 
that both equations possess periodic solutions for all positive values of \, and that these 
periodic solutions may be represented as power series in \, for small values of \, heuristic- 
ally by means of Lindstedt technique, or rigorously by means of the small-parameter 
Poincaré method, see Stoker, [5], or Lefschetz, [2]. 

These power series will not converge for all values of \. In the Van der Pol case, 
this is due to the fact that the equation does not possess a periodic solution for large 
negative values of \. In the case of the nonlinear spring, this is due to the existence of 
other singularities of the function. 

A question of some importance, from both the theoretical and computational point 
of view, is whether or not it is possible to obtain expansions which will be valid for all 
positive values of \, and, if not that, then at least valid for larger values of the positive 
parameter than are permitted by the power series expansion. 

We shall see that this is rigorously possible in the algebraic cases, and in the case 
of the linear integral equation, where the analytic character of the solution is readily 
obtained, and plausibly possible in the case of the Van der Pol and nonlinear spring 
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equations, where the analytic character of the solutions is less known. Actually, since 
the nonlinear spring equation may be solved explicitly in terms of elliptic functions, 
an exact investigation would not be too difficult in this case. However, we shall not 
include it here, since we are primarily interested in presenting and illustrating the basic 
principle in as simple and direct terms as possible. 

The fundamental idea is that of expanding the solution as a power series in another 
variable p which is itself a power series in A, p = @(A). 

We are essentially, then, presenting a method of analytic continuation. Many 
examples of this exist in the literature. The classical treatment of the hypergeometric 
equation exemplifies the general technique, and a particular version will be found in 
Shohat’s generally overlooked paper [4] on the Van der Pol equation. Unfortunately, 
the idea is a bit obscured by the introduction and subsequent elimination of a fictitious 
parameter. For a closely related discussion of summability properties of the Liouville- 
Neumann series obtained from the linear integral equation in (1.1d) see [1]. 

It is clear that it is easy to multiply examples of physical problems where this tech- 
nique may be of use. An example of particular importance, which we shall treat else- 
where, is the theory of shock waves of varying strength. As the strength increases, the 
ordinary perturbation techniques may be expected to be less and less effective. Here the 
physical background restricts us in a very natural way to the consideration of positive 
values. 

In passing, let us note that the Lighthill technique, see [3], may also be profitably 
modified in this manner in many cases. We shall discuss this topic in further detail 
at a subsequent time. 

2. Asimple algebraic example. Let us begin with an example which, if very simple, 
has the merit of illustrating the idea very clearly. Consider the problem of finding a 
power series development for the positive root of x” + x = \, where A > 0. Since we 
have an explicit formula for the solution in this case, we know that the power series 
for x in terms of \, z = \ — X* + --- , converges only for \ < 1/4. Let us perform the 
change of variable, p = A/(1 + 4A). Then 1 + 4A = 1/(1 — 4p) and 


z= —$4+ (1 — 4p)” = p+3p'4+---. (2.1) 


The radius of convergence in the p-plane is 1/4. Hence the series converges as long as 
| A/(1 + 4A) | < 1/4, which means that it converges for all positive i. 

The success of the method in this case is due to the fact that we know so much about 
the analytic character of x as a function of \. Suppose that we had merely set p equal 


to A/(1 + A). Then 


1 ea ae 


which means that the radius of convergence in the p-plane is 1/3. Hence as long as 
\/(1 +) | < 1/3 the expansion for z in terms of p will converge. Since 1/4/(1 + 1/4) = 


1/5 < 1/3, we see that the next expansion will allow larger positive values of \ than 


before. 
3. The solution of x + \Ax = c. Let us now turn to a discussion of the solution 


of the system in (1.1b). The following remarks will apply equally well to the solution 
of a linear integral equation, such as one in (1.1¢), and indeed are abstractly identical 


if we regard A as a positive operator. 
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The Liouville-Neumann solution of x + A\Ax = ¢, 
x=c—rA\Ac+NAC—::-, (3.1) 


obtained by iteration, converges within the circle, | \ | < 1/|Ay | where A, is the, or a, 
characteristic root of largest absolute value. 

If we take A to be a positive matrix, in the sense below*, the classical theorem of 
Perron asserts that there is a unique root of largest absolute value which is positive and 
simple. By a change of variable we may consider this root to be unity. 

To take advantage of the positivity of A, and the fact that at the moment we are 
interested only*in positive values of \, let us make the change of variable, p = /(1 + A), 
or \ = p/(1 — p). The equation for x becomes 


a(1 — p) + pAx =c— cp (3.2) 


with the power series solution 


z=c-— Acpt+ >. p(A — I)"(—Aeo). (3.3) 

Since the characteristic roots of A — J are X, — 1, k = 1, 2, --- , where i, are the 
roots of A, we see that the absolute value of any root of A — J is less than 2. Hence 
the radius of convergence of the series in (3.3), considered as a series in p is greater 
than 1/2. Since \/(1 + A) = 1/2 at A = 1, we see that in any case we have enlarged 


the set of positive values of \ which may be used for the series expansion in p. 

If, in addition to A being positive, it is also positive definite, then all of the charac- 
teristic roots of A will be positive, which means that 0 < 1 — , < 1, for all k. Hence 
the series in (3.3) will converge for | p | < 1, or! A/(1 + A) | < 1, which includes all 
positive values of X. 

In general, we would like to have the ratio of the root of largest absolute value to 
the root of next largest absolute value as large as possible, in order to speed up conver- 
gence. We can accomplish this by iterating in the usual fashion a finite number of times 
before introducing the above change of variable. For example, iterating twice, we derive 


x =c— rAc+d*A’ce — N*A*e (3.4) 


or 
z+ rA*x =c — \Ac + D°A’e, (3.5) 


where we now consider \*° to be the new parameter. 
4. The nonlinear spring. Let us now turn to the discussion of nonlinear differential 
equations, considering as a first example the equation 


a’ +a2+nd2* = 1, z(0) = 1, z’(0) = 0. (4.1) 


Following the ideas of Lindstedt, see [2], [4], [5], to obtain a series expansion for the 
periodic solution with the above normalized initial conditions, it is necessary to take 
account of the fact that nonlinearity affects.not only the amplitude, but also the fre- 
quency. This corresponds to the fact that while linear equations have fixed singularities, 
nonlinear equations in general have movable singularities, see [3]. 

Hence we make an initial change of variable, s = tv, where v is a parameter depending 


*A matrix A is called positive if ag; > 0, (i, 7 = 1, 2, -** , m). 
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upon \. The new equation is 


vac’ +2+d2° = 0, (4.2) 
which we write in the form 
(dv)*x’? + 7a + d®x* = O. (4.3) 
Let us set p = A/(1 + A), A = p/(1 — p), and 
vrA= p+ep +e, + ::: (4.4) 
since we want v to be 1 when A = O. We also set 
x = coss + pu,(s) + pur(s) + --- (4.5) 


and substitute in (4.1). The conditions that u, , uw, , and so on, be periodic functions of s 
with period 27 will determine the coefficients in the expansion in (4.4), as in the usual 
application of the Lindstedt method. The initial conditions, u;(0) = u{(0) = 0, will 
then determine the function u,(s) for 7 = 1, 2, 3, ---. 

It is useful to observe that the first approximations obtained from this method 
and the usual method will agree up to first degree terms in \. Consequently, a great 
deal of arithmetic work can be saved by using previous results. Furthermore, where 
previous perturbation techniques have carried the approximation to higher terms, 
these series can be transformed to yield the series above. 

The calculations indicated above are straightforward, and as always laborious. 
Computations very kindly performed by George Waters yield the following values 
for the coefficients in (4.4), 


1.37500 


1.66797 


Co 


C3 
Cc, = 1.91845 
c; = 2.14083 


(4.6) 


These coefficients yield a value of 1.24+ for »v when \ = 1, as compared with v = 1.28 
obtained from a REAC computation. 

Let us observe in passing that when the coefficients increase as gradually as those in 
(4.6), a crude extrapolation will yield values of c, and c; which will considerably decrease 
the error made in breaking off the series after five terms. Furthermore, if we are interested 
in the value \ = 1, it is probably better to extrapolate the values of c,(1/2)" based upon 
the known values n = 2, 3, 4, 5. 

5. The Van der Pol equation. Shohat expansion. Turning to the Van der Pol 
equation, the same change of variable as before yields the equation 


(Av)?ar’” + (Av)A*(a? — Ia’ + Ax = 0. (5.1) 
Setting 
X= p/(1—- pP)=ptet+et:::, 
vu = pt+ep +esp t+ep' +--:, (5.2) 
x = coss + u,(s) + u.(s) + --- , 
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the Lindstedt procedure furnishes the following values 
c¢=1, cs; = 15/16, c, = 13/16. 


From the following table, reproduced from Shohat’s paper cited above, it is tempting 
to conjecture that the Shohat series converges for all values of \ which are positive. 
If so, the series should be more widely known. 


r v computed using (5.2) v (Van der Pol) 
33 .98 .99 
1.0 .93 .90 
2.0 ae .78 
8.0 30 .39 
10.0 .30 31 


The Van der Pol values were obtained using graphical techniques. 
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ON MIDDLETON’S PAPER “SOME GENERAL RESULTS IN THE 
THEORY OF NOISE THROUGH NON-LINEAR DEVICES’’* 


By J. 8. SHIPMAN (Laboratory For Electronics, Inc., Boston, Mass.) 


As one of the central results of the title paper [1], Middleton obtained F,(t), the 
correlation function for the /th zone, as a function of the input correlation function ro 
in the case of the vth law half-wave rectification of narrow-band normal noise (see, e.g., 
his equations (7.14) and (7.15)). Unless one resorts to series evaluations, his formulas 
are not particularly suited for numerical computation as they stand, involving as they 
do hypergeometric functions which are not well tabulated. For purposes of calculation, 
then, a reduction of the hypergeometric functions occurring in the formulas to tabulated 
functions must ordinarily be effected, usually by applying the recursion relations among 
contiguous hypergeometric functions due to Gauss. 

When this reduction is accomplished, the hypergeometric functions in Middleton’s 
formulas are seen to be either polynomials in 7, or combinations of complete elliptic 
integrals of the first and second kind, provided » is an integer (see, e.g., Middleton’s 
equations (7.16) and (7.17)). These polynomials and combinations of elliptic integrals 
turn out to be, in every case so far examined, special cases of “Bennett functions” 
recently tabulated by the author and his colleagues [2, 3]. In the present note expressions 


*Received Sept. 24, 1954. Revised manuscript received Nov. 3, 1954. 
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for R,(t) in terms of Bennett functions for the important practical cases y = 1, 2 for 
several values of | and for 1 = 0 and v = 3, 4 are collected. The notation used is that of 
Middleton except for the symbol A&”? which is the Bennett function of the rth kind [3]. 


Thus 


Ro(t)ra1 = (8°p/2)[4A00 (r0)], (1) 
Ri(t)pa1 = (8’¥/2)(cos wet) Ady (ro), (2) 
R2{t),-1 = (8°y/2)(cos 2w.t) Aor (ro), (3) 
R(t). =0 (l=3,5,7-:°), (4) 
R,(t),-1 = (8°p/2)(cos 4w,t) Ads (To) ; (5) 
further 
Ro(t)r-2 = B¥'[}A00 (r0)], (6) 
R,(t),-2 = 8’ Y"(cos w,t) Asi (ro), (7) 
R.(t),-2 = B’ "(cos 2u,t)As2 (ro), (8) 
R;(t),-2 = B’Y"(cos 3w,t) Abs (To), (9) 
Ri(t)-.2.=0 (l= 4,6,8,---), (10) 
R;(t),-2 = B’W'(cos 5w,t) Abs (%); (11) 
and finally 
Ro(t),-3 = 38°Y*[3Aoo (r0)], (12) 
Ro(t) rma = 128°Y*[3A0co (ro) ]. (13) 


The functions A”; (r,) are tabulated directly for vy = 1, 2 [2, 3]; for »y > 3 recursion form- 
ulas are given which enable Bennett functions of higher kind to be obtained from those 
of lower kind. 

The author is indebted to E. Feuerstein and Dr. R. L. Sternberg for suggestions 
leading to the present note, and to Miss Jean Conway for assistance in the preparation 


of the manuscript. 
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MAXIMUM ACCELERATION IN TWO-DIMENSIONAL STEADY FLOWS 
OF AN IDEAL FLUID* 


By CHIA-SHUN YIH (Jowa Institute of Hydraulic Research) 


It is well known** that, in a domain free from singularities, the maximum speed in 
irrotational flows of an incompressible fluid occurs on the boundary of that domain. 
In this note, it will be proved that for two-dimensional irrotational flows of an incom- 
pressible fluid, the maximum magnitude of the acceleration in any singularity-free 
domain D must also occur on its boundary provided the flows are steady. 

The acceleration components in the z- and y-directions are, respectively, 


b = uu, + vu, , (1) 


c= w,+w,, (2) 
in which wu and v are the x- and y-components of the velocity, and the subscripts denote 
partial differentiation. The equation of continuity is 

uz +v, = 0 (3) 
and the equation expressing irrotationality is 

u, —v, = 0. (4) 
From Eqs. (1) to (4) it immediately follows that the square of the magnitude of the 
acceleration is given by 
a=) +c = g(ui+u), (5) 
in which g’ = wu? + v’. 

If @ and y are respectively the velocity potential and the stream function of the 
given flow, the complex potential for that flow is 
w(z) = o(2, y) + iW(z, y), 

in which w(z) is an analytic function free from singularities in the domain D. The well- 
known relation concerning the magnitude of the velocity is 


gq = w'(w’)* (6) 


in which the asterisk denotes the complex conjugate, and the primes denote ordinary 
differentiation. Since w(z) is analytic in D, its derivative 
| 


w'(z) = —u + w 


is also analytic in D and represents a potential flow with velocity components u, and wu, , 
concerning which the result analogous to (6) is 


2 vr sl) * (7) 
From (5) to (7), it follows that 


a = w'w'"(w'w"’)* = (w’/2)'[(w"/2)’]* = ai (8) 


*Received Nov. 1, 1954. This work is supported by the Office of Naval Research under Contract 
No. N8onr-500. 
**H. Lamb, Hydrodynamics, Dover, New York, 1945, p. 39 and p. 47. 
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in which gq, is the velocity magnitude in the irrotational flow corresponding to the new 
complex potential w’’/2. But it has been proved that the maximum value of gj in D 
occurs on its boundary, hence the same is true for a’. 

Since w’w” is analytic in D, the result for the acceleration could have been reached 
through 

a _ w'w''(w'w'’)* 
in (8)-by means of the well-known maximum-modulus theorem for analytic functions, 
without resorting to gq, . However, since the corresponding result for the speed is often 
reached without the use of the maximum-modulus theorem, the present proof may 
seem preferable to some. 

As a by-product of the proof, a means of evaluating the magnitude of the acceleration 
at any point in the field of flow is provided by (8), once w(z) is known. Furthermore, 
since the acceleration is proportional to the gradient of the dynamic pressure, it follows 
from the present result that the maximum dynamic-pressure gradient in D must occur 
on its boundary. 

The boundary of D is, of course, not necessarily a solid boundary. However, for 
ambiently uniform flows past arbitrarily shaped bodies, the maximum acceleration 
(or the maximum dynamic-pressure gradient) must occur on one of the solid boundaries. 
This can be seen by taking D as the entire domain outside the solid bodies. Its boundary 
then consists of the solid boundaries and a boundary at infinity. Since the acceleration 
at infinity is zero, the maximum acceleration must occur on one of the solid boundaries. 


AN ESTIMATE OF THE ERROR DUE TO THE TRUNCATED BOUNDARY 
IN THE NUMERICAL SOLUTION OF THE BLASIUS EQUATION 


By L. A. RUBEL! (UU. S. Naval Proving Ground*, Dahlgren, Va.) 


We seek a numerical solution of the Blasius differential equation: 
f''"(a) + fio f’'(x) = 0 
satisfying the split conditions 
0 = fQ) = f’'0, f'(e) = 2. 


The condition f’(*) = 2 is clearly unsuited for numerical procedures. In previous 
tabulations (e.g. [1], [2], [3]), the condition f’() = 2 is replaced by f’(M) = 2 where 
M is chosen so large that the second derivative becomes indistinguishable from zero, 
and further integration becomes impossible. That this procedure yields, in any sense, 
an approximation to the desired function has not, to the author’s knowledge, been 
established heretofore. 

This paper demonstrates that for large M, the numerical solution lies close to the 
true solution, gives a computable estimate of the deviation from the true solution, and 
provides an a priori criterion by which any desired accuracy may be achieved. 

1Received Nov. 12, 1954. The author is now at Cornell University. 

*The opinions and assertions expressed herein are the private ones of the writer and are not to be 
construed as official or reflecting the views of the Navy Department or the Naval Establishment at 


large. 
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Let f(x) satisfy 


f'’'+ff’ =0, f(0) = f’(O) 


I 
bo 


I 
—) 


f'(~) 
Let gu(x) = g(x) satisfy 
0, g'(M) = 2. 


II 


gy’ +99" =0, gO) = 9’) 


We observe some well known properties of f(x); namely that f(x) > 0 for x > 0; 
and therefore, that f’(z) increases steadily from 0 at x = 0 to2 atx = o. This may be 
demonstrated by observing that since f’(0) = O and f’(t) > O for some ¢ > 0, f’(z) 
must somewhere be an increasing function of z, and that f’’(z) must somewhere be 
positive. 

If f’’(x) is negative for some positive x, then f’’(x,) must be zero for some positive Zo - 
From the differential relation 


tf’ (Zo) + f (Xo) f’ (Xo) => 0 


we see that f’’(z)) = 0. Repeated differentiation of the differential equation demon- 
strates that f(z) = 0 for n > 2, and that, since f is analytic at x , f(x) is a poly- 
nomial, which is impossible. 

The function k’f’(kM) is zero for k = 0, increases with k, and is unbounded as k 
approaches ~. For some value of k, then, 


k*f"(kM) = 2. 


For this value of k, we see that g(x) = kf (kx) since kf (kx) satisfies the requirements 
that determine g(x) uniquely. It is clear that k > 1. 

Put g’'(@) = 2 + «(M) = 2k’ so that kh”? — 1 = 4¢(M). 

Let us now estimate | g(x) — f(x) | forz < M. 


| g(x) — f(x) | = | kf(ka) — f(x) | < | (& — If(kx) | + | (kx) — f(x) 
Since f(0) = O and f’(x) < 2 for z > 0, we have 
f(kx) < ke. 
By the law of the mean of differential calculus, 
flka) — f(x) = f'O(k — Dz for some &, x<st<ke. 
Since f’(é) < 2, 
f(kx) — f(x) | < 2k — Iz. 
Hence, 
| g(a) — f(x) | < 2k(k — lx + Ak — Va = Ak’ — 1)z. 
But k? — 1 = 4e(M), so that 
| g(x) — f(x) | < ae(M) < Me(M) for z< WM. 


We now estimate «(/). 


«(M) = g’(o) — g(M) = | g(x) de. 
M 
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Since 





or 
g’'(x) = g’"(M) exp | [ g(x) ac |, 


Now g(x) > g(M) for x > M, so that 


[ g(x) dx > g(M)(x — M), 
JM 


[ g’'(x) dx = g’"(M) / exp | - i) g(t) at] dx 
JM M M 


nn 
< g’’(M) [ exp [—g(M)(x — M)] dz = g(M) 
Hence 
y"M) 
(M) < Tip 


so that 


| f(x) — g(x) | < M ries for alla < M. 


Notice that this estimate involves only the tabulated function g and its derivatives 
so that the degree of approximation is directly computable. For example, the Emmons 
and Leigh tabulation [1] for M = 5 yields 


gi'(5) =0+4107,  g5(5) = 8.27922 + 107° 


so that 
| f(z) — g(x) |< $10, (x <5). 


It is possible to estimate beforehand how large M should be chosen to attain any 
desired accuracy, using the Emmons and Leigh tabulation. 
Let h(x) satisfy 


h’'"(x) + h(z)h’"(z) = 0, 
h(0) = h'(0) = 0, (5) = 2. 


Put h(x) = mf(mz). 
Let g(x) satisfy 
g’’"(x) + g(x)g’"(x) = 0, 
gO) = 90 =0, g(M)=2 (M>5). 


Put g(x) = kf (kx). Then 


ro = Eilts) 


Pr oer T gM) _ yy KKM /m) | 
| ote) — fa) |< Mn! = Ms rae ey 
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The tabulated values state: 
h’(1) = 1.25953 + 10°°, 


h’’(0) = 


— 


32823 + 10°°, 


h’(5) = 0+ 107°, 
h(5) = 8.27922 + 10°, 





lA 


4 10~°. 
8 


Now m* — 1 = 4¢(M) < 510°’ and m < 1+ 10°. Alsol < k < m. 
To estimate h(kM/m) and h’’(kM/m), use the crude estimation: 


h(x) > h’(A)(x — 1) 


and the representation 


h’’(x) = h’’(0) exp E [ h(t) a] < h’’(0) exp E [ h(t) ar], 
“ v1 


0 


After certain manipulations, which may be far from the best, the author has arrived 
at the following result: 
if M@ > 1+ 2(N + 1)”, then 


g(x) — f(x) | < 10~* tor zs. <= MM. 


Hence, M = 8 gives at least 10-place accuracy, and M = 10 gives at least 19-place 
accuracy. 

These last estimates are all crude. Once M is chosen, however, and the tabulation 
is made, Mg’’(M)/g(M) can be calculated easily, and measures the accuracy well. 


REFERENCES 
1. H. W. Emmons and D. Leigh, Tabulation of the Blasius function with blowing and suction, Harvard 
University, Division of Applied Sciences, Combustion Aerodynamics Laboratory, Interim 
Technical Report No. 9 
2. John C. Martin, On Magnus effects caused by the boundary layer displacement thickness on bodies 
of revolution at small angles of attack, Ballistic Research Laboratory, Report No. 870 
3. H. Schlichting, Boundary layer theory, Part I, Laminar Flows, NACA TM 1217, 1949 


ON THE EQUATIONS OF LINEARIZED CONICAL FLOW* 
By J. H. GIESE (Aberdeen Proving Ground, Maryland) 


In an irrotational flow field with velocity potential function (a, y, 2’) let u, v, w’ 
be velocity components at the point with rectangular coordinates z, y, 2’. If u, v, w’ are 
constant on rays through the origin the field is said to be conical [1, 2]. Then u, v, w’ 
are homogeneous of order zero in z, y, 2’, and ® is homogeneous of order one if an un- 
important arbitrary additive constant is neglected. Furthermore, if the flow field arises 
by slightly perturbing a uniform flow parallel to the 2’ axis at Mach number M, then 





*Received August 30, 1954. 
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® satisfies approximately 
0°b/dax*® + a°b/day? + (1 — M’) 0°’6/dz” = 0. 
The transformation z = i(M? — 1)7'”2’ takes this into 


a°b/dx? + a°b/dy’ + 3°b/dz? = 0. 


Thus the construction of linearized conical flow fields can be reduced to determination 
of homogeneous harmonic functions of order one. In this connection Poritsky [4, 5] has 
investigated homogeneous harmonic functions of order n, especially for n = 0 and 1, 
to obtain among other results well-known parametric representations of u, v, w’ for 
linearized conical flows in terms of analytic functions of a complex variable [2]. This 
note develops an alternative derivation of these formulas that is implicit in some of 
Busemann’s early work on conical flow [1]. 
Write (x, y, 2’) = zo(X, Y), where X = 2/z and Y = y/z. Then 


u = b/dx = d¢/AX, v = d8/dy = dg/aY, 
w = 00/02 = ¢ — uX —v!, (1) 
and Laplaces’s equation becomes 
(1 + X’) 0’y/aX? + 2XY d°o/aX AY + (1 + Y’) d’y/aY’ = 0. (2) 


It is easy to show that if u and v are functionally dependent zg must be linear in 2, y, z. 
Hence with no real loss of generality one may assume that u(X, Y) and v(X, Y) are 
functionally independent. Then exactly as in Busemann’s derivation of the equation 
for non-linearized conical flow [1], under the Legendre transformation (1), 


w, = dw/du = —X, w, = dw/dv = —Y, (3) 

and (2) becomes 
(1 + w?) 0°w/du? — 2w,w, d°w/du dv + (1 + w?) d’w/dv” = 0. (4) 
Now interpret u, v, w as rectangular coordinates in a “hodograph”’ space. By (4) the 
“‘hodograph”’ surface w = w(u, v) of a homogeneous harmonic function of order one is a 


minimal surface [1]. In accordance with the Weierstrass parametrization of minimal 
surfaces [3] 

i a 2 - +2 - 
u= Re- | (F° —G) dt, v=Re- 5? / (F° + G’) df, w = Re | FG dg (5) 


) 


for some analytic functions F and G of the complex variable ¢ = & + in. Conversely, 
for any analytic F(¢) and G(¢) for which u and v are functionally independent (5) defines 
a minimal surface. The so-called isothermal parameters £, » map the surface conformally 
onto the ¢-plane. 

To find the correspondence between the plane of Z = X, + 7Y and the ¢-plane note 
that (3) implies that the ray from (0, 0, 0) to (zx, y, z) is normal to the “hodograph”’ 
surface. Hence by (5) and the Cauchy-Riemann equations 

Z*F’ + 2FG — ZG’ = 0, 
where Z* = X — iY. This has the root 


Vie) = F/G = Z/[1 + (1 + 2Z2*%)"7) = (X + tY)+)/[1 + (1 + X? + OY’) #Y:~. 
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Then (1) and (5) imply 


-Re5[(V-V)af, =Re-Sif(Vt+V- df, w= Reso, 


o= 5K Isfer- v-yar—iv five vy arty], 


for some analytic f(¢). These comprise the desired results in the form exhibited by 
Poritsky, from which wu, v, w’ can easily be obtained as functions of z, y, 2’. It is customary 
to choose isothermal parameters such that V = ¢, so that in the supersonic case the 
interior of the Mach cone of the origin will map onto | ¢| < 1. 


REFERENCES 
A. Busemann, Aerodynamischer Auftrieb bei Ueberschallgeschwindigkeit, Luftfahrtforschung 12, 
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. H. Poritsky, Homogeneous harmonic functions, Quart. Appl. Math. 6, 379-388 (1949) 
. H. Poritsky, Linearized compressible flow, Quart. Appl. Math. 6, 389-404 (1949) 
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Introduction to aeronautical dynamics. By Manfred Rauscher. John Wiley & Sons, Inc., 
New York, and Chapman & Hall, Ltd., London. xv + 664 pp. $12.00. 


This book, as its unusual title indicates, has an unusual aim. It includes particle dynamics and 
rigid-body dynamics as well as the theory of fluid motion. The author obviously wishes to give a very 
detailed and thorough grounding in the fundamental ideas and notions in these fields. To this purpose, 
he provides detailed derivations of the basic laws; figures and graphs usually omitted in more concise 
presentations are generously included; basic examples are worked out in full detail; many original and 
stimulating examples and discussions are included, to help the beginning student overcome the usual 
basic difficulties. Thus, this book apparently has a higher aim than to be just a textbook; namely, to 
replace the teacher as well. Its features, described above, make it truly suitable for self-study, although 
there is no doubt that it will also be valuable in the classroom. 

As the various chapters proceed into more detailed and advanced matters, the author shows a 
certain preference for those special domains in which elegant solutions along classical lines can be ob- 
tained. For example, thin-airfoil theory is omitted in favor of an unusually complete presentation of 
airfoil generation by conformal mapping. Perhaps this is because otherwise the attempt to provide a 
book going into such detail would have led to a volume of unwieldy size and expense. At any rate, the 
result is a textbook of unusual character, reflecting, perhaps more than most books do, the special 


interests and predilections of its author. 
N. Rorr 


W. R. SEARS 


Einfiihrung in die Determinantentheorie einschliesslich der Fredhomschen Determinanten. 
By Gerhard Kowalewski. Walter de Gruyter & Co., Berlin, 1954. $7.16. 


The first edition of this well-known book was published in 1909. This fourth revised edition was 


prepared by the author, although publication was delayed due to his death in 1950. 
G. NEWELL 























